esei_EllipticCylinder.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluEllipticIntegralsEllipticCylinder.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_EllipticCylinder.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 eshelbySoluEllipticIntegralsEllipticCylinder::giveEllipticIntegrals(double J[13], double lambda, bool intpoint)
{
  // Elliptic cylinder: a1 > a2, a3 -> infinity
  if (I->a[0] <= I->a[1])
    {
      _warningg( "Values of semiaxes do not fit the condition for elliptic cylinder (a1 > a2, a3 is ignored). " );
      std::cout << "a1, a2, a3 = " << I->a[0] << ", " << I->a[1] << ", " << I->a[2];
    }

  double a1 = I->a[1];
  double a2 = I->a[2];
  // a3 is assumed to be infinite.

  CleanVector(J, 13);

  J[0] = 0.;

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

  // (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[9] = this->I23(a1, a2, J[2], J[3]); // I23
  J[8] = this->I22(a1, J[2], J[3], J[9], lambda); // I22

  J[10] = 0.0; // I31
  J[11] = J[9]; // I32
  J[12] = this->I33(a2, J[2], J[3], J[9], lambda); // I33

#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 eshelbySoluEllipticIntegralsEllipticCylinder::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] = (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+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
}

// I2 (corresponds to I1 in Mura's book, this cilinder is rotated and a1 is infinity)
double eshelbySoluEllipticIntegralsEllipticCylinder::I2(double a1, double a2, double lambda)
{
  return _4PI_ * a1 * a2 / (SQR(a2) - SQR(a1)) * (sqrt(SQR(a2) + lambda) / sqrt(SQR(a1) + lambda) - 1.0); // I2
}
// I3 (corresponds to I2 in Mura's book, this cilinder is rotated and a1 is infinity)
double eshelbySoluEllipticIntegralsEllipticCylinder::I3(double a1, double a2, double lambda)
{
  return _4PI_ * a1 * a2 / (SQR(a1) - SQR(a2)) * (sqrt(SQR(a1) + lambda) / sqrt(SQR(a2) + lambda) - 1.0); // I2
}
double eshelbySoluEllipticIntegralsEllipticCylinder::I22(double a1, double I2, double I3, double I23, double lambda)
{
  return ((I2 + I3) / (SQR(a1) + lambda) - I23) / 3.0;
}
double eshelbySoluEllipticIntegralsEllipticCylinder::I23(double a1, double a2, double I2, double I3)
{
  return (I3 - I2) / (SQR(a1) - SQR(a2));
}
double eshelbySoluEllipticIntegralsEllipticCylinder::I33(double a2, double I2, double I3, double I23, double lambda)
{
  return ((I2 + I3) / (SQR(a2) + lambda) - I23) / 3.0;
}

// Returns lambda for a given point coordinates
double eshelbySoluEllipticIntegralsEllipticCylinder::getLambda(const double a[3], double x1, double x2, double x3)
{
  // Solving quadratic equation.
  // a = 1.0
  double b = SQR(a[1]) + SQR(a[2]) - SQR(x2) - SQR(x3);
  double c = SQR(a[1]) * SQR(a[2]) - SQR(a[1]) * SQR(x3) - SQR(a[2]) * SQR(x2);
  double D = SQR(b) - 4.0 * c;
  if (D < 0.0) { _errorr("Discriminant is less than zero in eshelbySoluEllipticIntegralsEllipticCylinder::getLambda()"); }
  // Lambda is the largest positive root
  return (-b + sqrt(D)) / 2.0;
}


} // end of namespace mumech

/*end of file*/