esei_Cylinder.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluEllipticIntegralsCylinder.cpp
// author(s):      Jan Novak, Tomas Janda
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cmath>
#include <iostream>

#include "macros.h"
#include "esei_Cylinder.h"
#include "eshelbySoluLambda.h"
#include "legendreIntegrals.h"
#include "elementaryFunctions.h"
#include "inclusion.h"


namespace mumech {

// Function gives the values of Ferers-Dysons integrals
// Last edit:   14-Jan-2014
void eshelbySoluEllipticIntegralsCylinder::giveEllipticIntegrals(double J[13], double lambda, bool intpoint)
{
  // Cylinder: a2 = a3, a1 -> infinity
  if (I->a[1] != I->a[2])
    {
      _warningg( "Values of semiaxes do not fit the condition for circular cylinder (a2 = a3, a1 = infinity is ignored). " );
      std::cout << "a1, a2, a3 = " << I->a[0] << ", " << I->a[1] << ", " << I->a[2];
    }
  
  double a2 = I->a[1];
  CleanVector(J, 13);

  J[0] = 0.;

  J[1] = 0.0; // I1
  J[2] = this->I2(a2, lambda); // I2
  J[3] = J[2]; // I3 = I2

  // (Order changed due to sequential evaluation)
  J[4] = 0.0; // I11
  J[5] = 0.0; // I12
  J[6] = 0.0; // I13

  J[7] = 0.0; // I21
  J[8] = this->I22(a2, lambda); // I22
  J[9] = J[8]; // I23 = I22

  J[10] = 0.0; // I31
  J[11] = J[8]; // I32 = I22
  J[12] = J[8]; // I33 = I22

#ifdef TOM_DEBUG
  if (log)
    {
      std::cout << "J integrals" << std::endl;
      for (int i = 0; i < 13; i++)
    std::cout << "J[" << i << "] = " << J[i] << std::endl;
      std::cout << std::endl;
    }
#endif
}


// Function gives the derivatives of Ferers-Dysons integrals
// Last edit:   14-Jan-2014
// Flat shapes has zero drtivatives
void eshelbySoluEllipticIntegralsCylinder::giveDerivativesOfEllipticIntegrals(Point * point, bool intpoint)
{
  //if (point->position == PPF_INT_POINT) {
  //  CleanVector( point->dJi, 9 );
  //  CleanVector( point->dJij, 27 );
  //  CleanVector( point->ddJi, 27 );
  //  CleanVector( point->ddJij, 81 );
  //}
  //else if (point->position == PPF_EXT_POINT) {
  // NUMERICAL DERIVATIVES (using REAL ARGUMENTS)

  double ndiff = I->ndiff_1;

  // Get perturbated points and their lambdas
  double lambdas[NUM_PERTURB];
  this->getPerturbatedLambdas(lambdas, point->loc_x);

  // Create and fill a two dimensional array containing values of J-integrals at perturbated points.
  const int NUM_ELLIP_INT = 13;
  double J_pert[NUM_PERTURB][NUM_ELLIP_INT]; // All perturbations of all elliptical integrals.
      
  this->log = false;
  for (int i = 0; i < NUM_PERTURB; i++) // For each perturbated point
    this->giveEllipticIntegrals(J_pert[i], lambdas[i], intpoint); // Get elliptical integrals for lambda of perturbated point.
  this->log = true;

  // Compute numerial derivatives of first and second order of J_i integrals.
  for (int i = 0; i < 3; i++) // For each elliptic integral Ji
    {
      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Ji_j = dJi[3*i+j-4], i,j=1..3.
      //point->dJi[3*i] =   (J_pert[_Xp_][i+1] - J_pert[_Xm_][i+1]) / (2.0 * ndiff);
      point->dJi[3*i+1] = (J_pert[_Yp_][i+1] - J_pert[_Ym_][i+1]) / (2.0 * ndiff);
      point->dJi[3*i+2] = (J_pert[_Zp_][i+1] - J_pert[_Zm_][i+1]) / (2.0 * ndiff);

      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Ji_jk = ddJi[9*i+3*j+k-13], i,j,k=1..3.
      //point->ddJi[9*i] =   (J_pert[_Xp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Xm_][i+1]) / (ndiff * ndiff); // d^2f(x,y,z)/dx^2 = (f(x+h,y,z) - 2*f(x,y,z) + f(x-h,y,z)) / h^2
      //point->ddJi[9*i+1] = (J_pert[_XpYp_][i+1] - J_pert[_XpYm_][i+1] - J_pert[_XmYp_][i+1] + J_pert[_XmYm_][i+1]) / (4.0*ndiff*ndiff); // d^2f(x,y,z)/dxdy = (f(x+h,y+h,z) - f(x+h,y-h,z) - f(x-h,y+h,z) + f(x-h,y-h,z)) / (4*h^2)
      //point->ddJi[9*i+2] = (J_pert[_XpZp_][i+1] - J_pert[_XpZm_][i+1] - J_pert[_XmZp_][i+1] + J_pert[_XmZm_][i+1]) / (4.0*ndiff*ndiff);

      //point->ddJi[9*i+3] = point->ddJi[9*i+1];
      point->ddJi[9*i+4] = (J_pert[_Yp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Ym_][i+1]) / (ndiff * ndiff);
      point->ddJi[9*i+5] = (J_pert[_YpZp_][i+1] - J_pert[_YpZm_][i+1] - J_pert[_YmZp_][i+1] + J_pert[_YmZm_][i+1]) / (4.0*ndiff*ndiff);

      //point->ddJi[9*i+6] = point->ddJi[9*i+2];
      point->ddJi[9*i+7] = point->ddJi[9*i+5];
      point->ddJi[9*i+8] = (J_pert[_Zp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Zm_][i+1]) / (ndiff * ndiff);
    }

  // Compute numerial derivatives of first and second order of J_ij integrals.
  for (int i = 0; i < 9; i++) // For each elliptic integral Jij
    {
      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Jij_k = dJij[9*k+3*i+j-13], i,j,k=1..3.
      //point->dJij[i] =    (J_pert[_Xp_][i+4] - J_pert[_Xm_][i+4]) / (2.0 * ndiff);
      point->dJij[i+9] = (J_pert[_Yp_][i+4] - J_pert[_Ym_][i+4]) / (2.0 * ndiff);
      point->dJij[i+18] = (J_pert[_Zp_][i+4] - J_pert[_Zm_][i+4]) / (2.0 * ndiff);

      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Jij_kl = ddJi[27*i+9*j+3*k+l-40] i,j,k,l=1..3.
      //point->ddJij[9*i] =   (J_pert[_Xp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Xm_][i+4]) / (ndiff * ndiff);
      //point->ddJij[9*i+1] = (J_pert[_XpYp_][i+4] - J_pert[_XpYm_][i+4] - J_pert[_XmYp_][i+4] + J_pert[_XmYm_][i+4]) / (4.0*ndiff*ndiff);
      //point->ddJij[9*i+2] = (J_pert[_XpZp_][i+4] - J_pert[_XpZm_][i+4] - J_pert[_XmZp_][i+4] + J_pert[_XmZm_][i+4]) / (4.0*ndiff*ndiff);

      //point->ddJij[9*i+3] = point->ddJij[9*i+1];
      point->ddJij[9*i+4] = (J_pert[_Yp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Ym_][i+4]) / (ndiff * ndiff);
      point->ddJij[9*i+5] = (J_pert[_YpZp_][i+4] - J_pert[_YpZm_][i+4] - J_pert[_YmZp_][i+4] + J_pert[_YmZm_][i+4]) / (4.0*ndiff*ndiff);

      //point->ddJij[9*i+6] = point->ddJij[9*i+2];
      point->ddJij[9*i+7] = point->ddJij[9*i+5];
      point->ddJij[9*i+8] = (J_pert[_Zp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Zm_][i+4]) / (ndiff * ndiff);
    }
  //}
  //else _errorr("point position is not defined");
    
#ifdef TOM_DEBUG
  std::cout << "dJi integrals" << std::endl;
  for (int i = 0; i < 9; i++)
    std::cout << "dJi[" << i << "] = " << point->dJi[i] << std::endl;
  std::cout << std::endl;
  std::cout << "dJij integrals" << std::endl;
  for (int i = 0; i < 27; i++)
    std::cout << "dJij[" << i << "] = " << point->dJij[i] << std::endl;
  std::cout << std::endl;
  std::cout << "ddJi integrals" << std::endl;
  for (int i = 0; i < 27; i++)
    std::cout << "ddJi[" << i << "] = " << point->ddJi[i] << std::endl;
  std::cout << std::endl;
  std::cout << "ddJij integrals" << std::endl;
  for (int i = 0; i < 81; i++)
    std::cout << "ddJij[" << i << "] = " << point->ddJij[i] << std::endl;
  std::cout << std::endl;
#endif
}

double eshelbySoluEllipticIntegralsCylinder::I2(double a2, double lambda)
{
  return 2.0 * PI * SQR(a2) / (SQR(a2) + lambda); // I2
}

double eshelbySoluEllipticIntegralsCylinder::I22(double a2, double lambda)
{
  return PI * SQR(a2) / (SQR(a2) + lambda) / (SQR(a2) + lambda) ;
}

// Returns lambda for a given point coordinates
double eshelbySoluEllipticIntegralsCylinder :: getLambda (const double a[3], double x1, double x2, double x3)
{
  // The cubic equation simplify to this expression for circular cylinder.
  return SQR(x2) + SQR(x3) - SQR(a[1]);
}


} // end of namespace mumech

/*end of file*/