AnalyticalFun.hpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM 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 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM 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 MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
namespace FractureMechanics {
 
struct AnalyticalOptions {
  double aLpha;
  bool optionsInitialised;
  PetscErrorCode getOptions() {
 
    MoFEMFunctionBeginHot;
    ierr =
        PetscOptionsBegin(PETSC_COMM_WORLD, "",
                          "Get analytical boundary conditions options", "none");
    CHKERRQ(ierr);
    ierr = PetscOptionsScalar("-analytical_alpha", "crack angle", "", aLpha,
                              &aLpha, PETSC_NULL);
    CHKERRQ(ierr);
    PetscBool flg;
    double fraction = 1;
    ierr = PetscOptionsScalar(
        "-analytical_alpha_fraction_pi",
        "crack angle given as fraction of Pi aLpha = M_PI/fraction", "",
        fraction, &fraction, &flg);
    CHKERRQ(ierr);
    if (flg == PETSC_TRUE) {
      aLpha = M_PI / fraction;
    }
    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);
    optionsInitialised = true;
    MoFEMFunctionReturnHot(0);
  }
  AnalyticalOptions() : aLpha(0), optionsInitialised(false) {}
};
 
#ifndef __ANALITICAL_TRACTION__
#define __ANALITICAL_TRACTION__
 
struct AnalyticalForces : public NeumannForcesSurface::MethodForAnalyticalForce,
                          public AnalyticalOptions {
 
  PetscErrorCode getForce(const EntityHandle ent, const VectorDouble3 &coords,
                          const VectorDouble3 &normal, VectorDouble3 &forces) {
 
    MoFEMFunctionBeginHot;
 
    // initailise options
    if (!optionsInitialised) {
      ierr = getOptions();
      CHKERRQ(ierr);
    }
 
    const double x = coords[0];
    const double y = coords[1];
    const double z = coords[2];
    double &fx = forces[0];
    double &fy = forces[1];
    double &fz = forces[2];
 
    const double n0 = normal[0] / sqrt(pow(normal[0], 2) + pow(normal[1], 2) +
                                       pow(normal[2], 2));
    const double n1 = normal[1] / sqrt(pow(normal[0], 2) + pow(normal[1], 2) +
                                       pow(normal[2], 2));
    const double n2 = normal[2] / sqrt(pow(normal[0], 2) + pow(normal[1], 2) +
                                       pow(normal[2], 2));
 
    const double pi = M_PI;
    const double alpha = aLpha;
    const double theta = atan2(z, x * sin(alpha) + y * cos(alpha));
 
    const double rxy =
        sqrt(pow(x * sin(alpha) + y * cos(alpha), 2) + pow(z, 2));
    const double sref1 = 1 / sqrt(2 * pi * rxy) * cos(theta / 2.) *
                         (1. - sin(theta / 2.) * sin(3. * theta / 2.));
    const double sref2 = 1 / sqrt(2 * pi * rxy) * cos(theta / 2.) *
                         (1. + sin(theta / 2.) * sin(3. * theta / 2.));
    const double sref3 = 1 / sqrt(2 * pi * rxy) * sin(theta / 2.) *
                         cos(theta / 2.) * cos(3. * theta / 2.);
 
    const double scal1 = sref1 * pow(sin(alpha), 2);
    const double scal2 = sref1 * pow(cos(alpha), 2);
    const double scal3 = sref2;
    const double scal4 = sref1 * sin(alpha) * cos(alpha);
    const double scal5 = sref3 * sin(alpha);
    const double scal6 = sref3 * cos(alpha);
 
    // This is just example
    fx = scal1 * n0 + scal4 * n1 + scal5 * n2;
    fy = scal4 * n0 + scal2 * n1 + scal6 * n2;
    fz = scal5 * n0 + scal6 * n1 + scal3 * n2;
 
    fx *= -1;
    fy *= -1;
    fz *= -1;
 
    MoFEMFunctionReturnHot(0);
  }
};
 
#endif //__ANALITICAL_TRACTION__
 
#ifndef __ANALITICAL_DISPLACEMENT__
#define __ANALITICAL_DISPLACEMENT__
 
struct AnalyticalDisp : public AnalyticalOptions {
 
  vector<VectorDouble> dIsp;
 
  virtual vector<VectorDouble> &operator()(const double X, const double Y,
                                           const double Z) {
 
    // initailise options
    if (!optionsInitialised) {
      getOptions();
    }
 
    dIsp.resize(1);
    dIsp[0].resize(3);
 
    double &x = dIsp[0][0];
    double &y = dIsp[0][1];
    double &z = dIsp[0][2];
 
    // Current positions is given by
 
    const double E = 1.0e5;
    const double nu = 0;
    const double pi = M_PI;
    const double alpha = aLpha;
 
    const double theta = atan2(Z, X * sin(alpha) + Y * cos(alpha));
 
    const double ux =
        (1 + nu) / E *
        sqrt(sqrt(pow(X * sin(alpha) + Y * cos(alpha), 2) + pow(Z, 2)) /
             (2 * pi)) *
        cos(theta / 2.) * (3 - 4 * nu - cos(theta)) *
        sin(alpha); // simply ux = f(Z) = Z*1e-3
    const double uy =
        (1 + nu) / E *
        sqrt(sqrt(pow(X * sin(alpha) + Y * cos(alpha), 2) + pow(Z, 2)) /
             (2 * pi)) *
        cos(theta / 2.) * (3 - 4 * nu - cos(theta)) * cos(alpha);
    const double uz =
        (1 + nu) / E *
        sqrt(sqrt(pow(X * sin(alpha) + Y * cos(alpha), 2) + pow(Z, 2)) /
             (2 * pi)) *
        sin(theta / 2.) * (3 - 4 * nu - cos(theta));
 
    x = X + ux;
    y = Y + uy;
    z = Z + uz;
 
    return dIsp;
  }
 
};
 
// struct AnalyticalForces: public NeumannForcesSurface::MethodForAnalyticalForce {
//
//   PetscErrorCode getForce(
//     const EntityHandle ent,
//     const VectorDouble3 &coords,
//     const VectorDouble3 &normal,
//     VectorDouble3 &forces
//   ) {
//     MoFEMFunctionBeginHot;
//
//     const double x = coords[0];
//     const double y = coords[1];
//     const double z = coords[2];
//     double &fx = forces[0];
//     double &fy = forces[1];
//     double &fz = forces[2];
//
//     // This is just example
//     fx = z*1000;
//     fy = 0;
//     fz = 0;
//
//     // cerr << forces << endl;
//
//     MoFEMFunctionReturnHot(0);
//   }
//
// };
//
// #endif //__ANALITICAL_TRACTION__
//
// #ifndef __ANALITICAL_DISPLACEMENT__
// #define __ANALITICAL_DISPLACEMENT__
//
// struct AnalyticalDisp {
//
//   vector<VectorDouble > dIsp;
//
//   virtual vector<VectorDouble >& operator()(
//     const double X, const double Y, const double Z
//   ) {
//
//     dIsp.resize(1);
//     dIsp[0].resize(3);
//
//     double &x = dIsp[0][0];
//     double &y = dIsp[0][1];
//     double &z = dIsp[0][2];
//
//     // Current positions is given by
//
//     const double ux = Z*1e-3; // simply ux = f(Z) = Z*1e-3
//     const double uy = 0;
//     const double uz = 0;
//
//     x = X+ux;
//     y = Y+uy;
//     z = Z+uz;
//
//     return dIsp;
//
//   }
//
// };
 
}
 
#endif //__ANALITICAL_DISPLACEMENT__