v0.10.0
Remodeling.cpp

Integration points are taken from OOFEM http://www.oofem.org/resources/doc/oofemrefman/microplanematerial_8C_source.html

/**
* \file Remodeling.cpp
* \example Remodeling.cpp
*
* \brief Integration points are taken from OOFEM
* <http://www.oofem.org/resources/doc/oofemrefman/microplanematerial_8C_source.html>
*
*/
/*
* 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/>. */
#include <boost/program_options.hpp>
using namespace std;
namespace po = boost::program_options;
#include <BasicFiniteElements.hpp>
using namespace MoFEM;
#include <Remodeling.hpp>
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
#define TENSOR1_VEC_PTR(VEC) &VEC[0], &VEC[1], &VEC[2]
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT) \
&MAT(0, 0), &MAT(0, 1), &MAT(0, 2), &MAT(0, 3), &MAT(0, 4), &MAT(0, 5), \
&MAT(1, 0), &MAT(1, 1), &MAT(1, 2), &MAT(1, 3), &MAT(1, 4), &MAT(1, 5), \
&MAT(2, 0), &MAT(2, 1), &MAT(2, 2), &MAT(2, 3), &MAT(2, 4), &MAT(2, 5), \
&MAT(3, 0), &MAT(3, 1), &MAT(3, 2), &MAT(3, 3), &MAT(3, 4), &MAT(3, 5), \
&MAT(4, 0), &MAT(4, 1), &MAT(4, 2), &MAT(4, 3), &MAT(4, 4), &MAT(4, 5), \
&MAT(5, 0), &MAT(5, 1), &MAT(5, 2), &MAT(5, 3), &MAT(5, 4), &MAT(5, 5)
#define TENSOR4_MAT_PTR(MAT) &MAT(0, 0), MAT.size2()
#define TENSOR2_MAT_PTR(MAT) \
&MAT(0, 0), &MAT(1, 0), &MAT(2, 0), &MAT(3, 0), &MAT(4, 0), &MAT(5, 0), \
&MAT(6, 0), &MAT(7, 0), &MAT(8, 0)
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT) \
&MAT(0, 0), &MAT(0, 1), &MAT(0, 2), &MAT(0, 3), &MAT(0, 4), &MAT(0, 5)
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC) \
&VEC[0], &VEC[1], &VEC[2], &VEC[3], &VEC[4], &VEC[5]
namespace BoneRemodeling {
MoFEMErrorCode Remodeling::CommonData::getParameters() {
PetscBool flg_config = PETSC_FALSE;
is_atom_testing = PETSC_FALSE;
with_adol_c = PETSC_FALSE;
char my_config_file_name[255];
double young_modulus = 5.0e+9; // Set here some typical value for bone 5 GPa
double poisson_ratio = 0.3; // Set here some typical value for bone
c = 4.0e-5; // Set here some typical value for bone. days/ m2 // this
// parameters governs how fast we can achieve equilibrium
m = 3.25; // Set some parameter typical for bone
n = 2.25; // Set Some parameter typical for bone
rHo_ref = 1000; // Set Some parameter typical for bone. kg/m3
pSi_ref = 2.75e4; // Set Some parameter typical for bone. energy / volume unit
// J/m3 = Pa
rHo_max = rHo_ref + 0.5 * rHo_ref;
rHo_min = rHo_ref - 0.5 * rHo_ref;
b = 0;
R0 = 0.0; // Set Some parameter typical for bone.
cUrrent_psi = 0.0;
cUrrent_mass = 0.0;
nOremodellingBlock = false;
less_post_proc = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling parameters",
"none");
// config file
CHKERR PetscOptionsString("-my_config", "configuration file name", "",
"my_config.in", my_config_file_name, 255,
&flg_config);
if (flg_config) {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Config file: %s loaded.\n",
my_config_file_name);
try {
ifstream ini_file(my_config_file_name);
if (!ini_file.is_open()) {
SETERRQ(PETSC_COMM_SELF, 1, "*** -my_config does not exist ***");
}
po::variables_map vm;
po::options_description config_file_options;
// std::cout << my_config_file_name << std::endl;
config_file_options.add_options()(
"young_modulus", po::value<double>(&young_modulus)->default_value(1));
config_file_options.add_options()(
"poisson_ratio", po::value<double>(&poisson_ratio)->default_value(1));
config_file_options.add_options()(
"c", po::value<double>(&c)->default_value(1));
config_file_options.add_options()(
"m", po::value<double>(&m)->default_value(1));
config_file_options.add_options()(
"n", po::value<double>(&n)->default_value(1));
config_file_options.add_options()(
"rHo_ref", po::value<double>(&rHo_ref)->default_value(1));
config_file_options.add_options()(
"pSi_ref", po::value<double>(&pSi_ref)->default_value(1));
config_file_options.add_options()(
"R0", po::value<double>(&R0)->default_value(1));
po::parsed_options parsed =
parse_config_file(ini_file, config_file_options, true);
store(parsed, vm);
po::notify(vm);
} catch (const std::exception &ex) {
std::ostringstream ss;
ss << ex.what() << std::endl;
SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
}
}
CHKERR PetscOptionsScalar("-young_modulus", "get young modulus", "",
young_modulus, &young_modulus, PETSC_NULL);
CHKERR PetscOptionsScalar("-poisson_ratio", "get poisson_ratio", "",
poisson_ratio, &poisson_ratio, PETSC_NULL);
CHKERR PetscOptionsScalar("-c", "density evolution (growth) velocity [d/m^2]",
"", c, &c, PETSC_NULL);
CHKERR PetscOptionsScalar("-m", "algorithmic exponent", "", m, &m,
PETSC_NULL);
CHKERR PetscOptionsScalar("-n", "porosity exponent", "", n, &n, PETSC_NULL);
CHKERR PetscOptionsInt("-b", "bell function exponent", "", b, &b, PETSC_NULL);
CHKERR PetscOptionsScalar("-rho_ref", "reference bone density", "", rHo_ref,
&rHo_ref, PETSC_NULL);
CHKERR PetscOptionsScalar("-rho_max", "reference bone density", "", rHo_max,
&rHo_max, PETSC_NULL);
CHKERR PetscOptionsScalar("-rho_min", "reference bone density", "", rHo_min,
&rHo_min, PETSC_NULL);
CHKERR PetscOptionsScalar("-psi_ref", "reference energy density", "", pSi_ref,
&pSi_ref, PETSC_NULL);
CHKERR PetscOptionsScalar("-r0", "mass source", "", R0, &R0, PETSC_NULL);
CHKERR PetscOptionsBool("-my_is_atom_test",
"is used with testing, exit with error when diverged",
"", is_atom_testing, &is_atom_testing, PETSC_NULL);
CHKERR PetscOptionsBool("-less_post_proc",
"is used to reduce output file size", "",
less_post_proc, &less_post_proc, PETSC_NULL);
CHKERR PetscOptionsBool(
"-with_adolc",
"calculate the material tangent with automatic differentiation", "",
with_adol_c, &with_adol_c, PETSC_NULL);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Young's modulus E[Pa]: %4.3g\n",
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Poisson ratio nu[-] %4.3g\n",
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Lame coefficient lambda[Pa]: %4.3g\n", lambda);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Lame coefficient mu[Pa]: %4.3g\n", mu);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Density growth velocity [d/m2] %4.3g\n", c);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Algorithmic exponent m[-]: %4.3g\n", m);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Porosity exponent n[-]: %4.3g\n", n);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Reference density rHo_ref[kg/m3]: %4.3g\n", rHo_ref);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Reference energy pSi_ref[Pa]: %4.3g\n", pSi_ref);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Mass conduction R0[kg/m3s]: %4.3g\n", R0);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Bell function coefficent b[-]: %4.3g\n", b);
if (b != 0) {
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Max density rHo_max[kg/m3]: %4.3g\n", rHo_max);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Min density rHo_min[kg/m3]: %4.3g\n", rHo_min);
}
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
}
inline double Remodeling::CommonData::getCFromDensity(const double &rho) {
double mid_rho = 0.5 * (rHo_max + rHo_min);
double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
return 1 / (1 + (b != 0) * pow(var_h, 2 * b));
}
inline double Remodeling::CommonData::getCFromDensityDiff(const double &rho) {
double mid_rho = 0.5 * (rHo_max + rHo_min);
double var_h = (rho - mid_rho) / (rHo_max - mid_rho);
return (b != 0) * (-2) * b * pow(var_h, 2 * b - 1) /
((rHo_max - mid_rho) * pow(pow(var_h, 2 * b) + 1, 2));
}
OpGetRhoTimeDirevative::OpGetRhoTimeDirevative(
Remodeling::CommonData &common_data)
"RHO", UserDataOperator::OPROW),
commonData(common_data) {}
OpGetRhoTimeDirevative::doWork(int side, EntityType type,
const int nb_gauss_pts = data.getN().size1();
if (!nb_gauss_pts)
const int nb_base_functions = data.getN().size2();
auto base_function = data.getFTensor0N();
commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
if (type == MBVERTEX) {
commonData.data.vecRhoDt.clear();
}
int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
double *array;
CHKERR VecGetArray(getFEMethod()->ts_u_t, &array);
auto drho_dt = getFTensor0FromVec(commonData.data.vecRhoDt);
for (int gg = 0; gg < nb_gauss_pts; gg++) {
int bb = 0;
for (; bb < nb_dofs; bb++) {
const double field_data = array[data.getLocalIndices()[bb]];
drho_dt += field_data * base_function;
++base_function;
}
// It is possible to have more base functions than dofs
for (; bb != nb_base_functions; bb++)
++base_function;
++drho_dt;
}
CHKERR VecRestoreArray(getFEMethod()->ts_u_t, &array);
}
OpCalculateStress::OpCalculateStress(Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {
// Set that this operator is executed only for vertices on element
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
OpCalculateStress::doWork(int side, EntityType type,
if (type != MBVERTEX)
if (!commonData.data.gradDispPtr) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Gradient at integration points of displacement is not calculated");
}
auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
if (!commonData.data.rhoPtr) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Density at integration points is not calculated");
}
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
commonData.data.vecPsi.resize(nb_gauss_pts, false);
auto psi = getFTensor0FromVec(commonData.data.vecPsi);
commonData.data.vecR.resize(nb_gauss_pts, false);
auto R = getFTensor0FromVec(commonData.data.vecR);
commonData.data.matInvF.resize(9, nb_gauss_pts, false);
auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
commonData.data.vecDetF.resize(nb_gauss_pts, false);
auto det_f = getFTensor0FromVec(commonData.data.vecDetF);
MatrixDouble &matS = commonData.data.matS;
matS.resize(nb_gauss_pts, 6, false);
commonData.data.matP.resize(9, nb_gauss_pts, false);
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
MatrixDouble &matF = commonData.data.matGradientOfDeformation;
matF.resize(9, nb_gauss_pts, false);
const double c = commonData.c;
const double n = commonData.n;
const double m = commonData.m;
const double rho_ref = commonData.rHo_ref;
const double psi_ref = commonData.pSi_ref;
const double lambda = commonData.lambda;
const double mu = commonData.mu;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get deformation gradient
F(i, I) = grad(i, I);
F(0, 0) += 1;
F(1, 1) += 1;
F(2, 2) += 1;
CHKERR invertTensor3by3(F, det_f, invF);
const double log_det = log(det_f);
psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
psi *= 0.5 * mu;
psi += 0.5 * lambda * log_det * log_det;
const double coef = lambda * log_det - mu;
P(i, J) = mu * F(i, J) + coef * invF(J, i);
S(i, I) = invF(i, J) ^ P(J, I);
const double kp = commonData.getCFromDensity(rho);
R = c * kp * (pow(rho / rho_ref, n - m) * psi - psi_ref);
++det_f;
++invF;
++grad;
++rho;
++psi;
++S;
++P;
++R;
++F;
}
}
OpCalculateStressTangentWithAdolc::OpCalculateStressTangentWithAdolc(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {
// Set that this operator is executed only for vertices on element
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
MoFEMErrorCode OpCalculateStressTangentWithAdolc::doWork(
int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
if (type != MBVERTEX)
auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
vecC.resize(6, false);
MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
dS_dC.resize(6, 6 * nb_gauss_pts, false);
dS_dC.clear();
MatrixDouble &matS = commonData.data.matS;
matS.resize(nb_gauss_pts, 6, false);
MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
dP_dF.resize(81, nb_gauss_pts, false);
MatrixDouble &matF = commonData.data.matGradientOfDeformation;
matF.resize(9, nb_gauss_pts, false);
commonData.data.matP.resize(9, nb_gauss_pts, false);
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
commonData.data.vecPsi.resize(nb_gauss_pts, false);
auto psi = getFTensor0FromVec(commonData.data.vecPsi);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get deformation gradient
F(i, I) = grad(i, I);
F(0, 0) += 1;
F(1, 1) += 1;
F(2, 2) += 1;
// Calculate Green defamation Tensor0
// ^ that means multiplication of Tensor 2 which result in symmetric Tensor
// rank 2.
C(I, J) = F(i, I) ^ F(i, J);
// Calculate free energy
CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
commonData, C, psi);
int r_g = ::gradient(commonData.tAg, 6, &vecC[0], &matS(gg, 0));
if (r_g < 0) {
// That means that energy function is not smooth and derivative
// can not be calculated,
SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"ADOL-C function evaluation with error");
}
// Scale energy density and 2nd Piola Stress
S(0, 0) *= 2;
S(1, 1) *= 2;
S(2, 2) *= 2;
// Calculate 1st Piola-Stress tensor
P(i, J) = F(i, I) * S(I, J);
const int is = 6 * gg;
double *hessian_ptr[6] = {&dS_dC(0, is), &dS_dC(1, is), &dS_dC(2, is),
&dS_dC(3, is), &dS_dC(4, is), &dS_dC(5, is)};
int r_h = ::hessian(commonData.tAg, 6, &vecC[0], hessian_ptr);
if (r_h < 0) {
// That means that energy function is not smooth and derivative
// can not be calculated,
SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
"ADOL-C function evaluation with error");
}
// FIXME: TODO: Here tensor 4 with two minor symmetries is used, since
// \psi_{,ijkl} = \psi_{,klij} = \psi_{,lkji} = ... there is additional
// major symmetry. That is 21 independent components not 36 how is now
// implemented using Tenso4_ddg.
//
// That way we have to fill off-diagonal part for dS_dC. In future that
// need to be fixed by implementation of Tenso4 with additional major
// symmetries. Pleas consider doing such implementation by yourself. You
// have support of MoFEM developing team to guide you through this process.
//
// This deteriorate efficiency.
//
// TODO: Implement matrix with more symmetries.
//
for (int ii = 0; ii < 6; ii++) {
for (int jj = 0; jj < ii; jj++) {
dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
}
}
// As result that symmetry of deformation gradient C is used, terms which
// have mixed indices are two time to big, and terms which has two times
// midxed index are four times to big. Those terms are not multiplied.
D1(0, 0, 0, 0) *= 4;
D1(1, 1, 1, 1) *= 4;
D1(2, 2, 2, 2) *= 4;
D1(0, 0, 1, 1) *= 4;
D1(1, 1, 0, 0) *= 4;
D1(1, 1, 2, 2) *= 4;
D1(2, 2, 1, 1) *= 4;
D1(0, 0, 2, 2) *= 4;
D1(2, 2, 0, 0) *= 4;
D1(0, 0, 0, 1) *= 2;
D1(0, 0, 0, 2) *= 2;
D1(0, 0, 1, 2) *= 2;
D1(0, 1, 0, 0) *= 2;
D1(0, 2, 0, 0) *= 2;
D1(1, 2, 0, 0) *= 2;
D1(1, 1, 0, 1) *= 2;
D1(1, 1, 0, 2) *= 2;
D1(1, 1, 1, 2) *= 2;
D1(0, 1, 1, 1) *= 2;
D1(0, 2, 1, 1) *= 2;
D1(1, 2, 1, 1) *= 2;
D1(2, 2, 0, 1) *= 2;
D1(2, 2, 0, 2) *= 2;
D1(2, 2, 1, 2) *= 2;
D1(0, 1, 2, 2) *= 2;
D1(0, 2, 2, 2) *= 2;
D1(1, 2, 2, 2) *= 2;
D2(i, J, k, L) = (F(i, I) * (D1(I, J, K, L) * F(k, K)));
// D2(i, J, k, L) += I(i, k) * S(J, L);
++grad;
++S;
++F;
++D1;
++D2;
++P;
++psi;
}
}
std::vector<EntityHandle> &map_gauss_pts,
Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data), postProcMesh(post_proc_mesh),
mapGaussPts(map_gauss_pts) {
// Set that this operator is executed only for vertices on element
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
OpPostProcStress::doWork(int side, EntityType type,
if (type != MBVERTEX)
double def_VAL[9];
bzero(def_VAL, 9 * sizeof(double));
Tag th_piola2;
CHKERR postProcMesh.tag_get_handle("Piola_Stress2", 9, MB_TYPE_DOUBLE,
th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
Tag th_psi;
CHKERR postProcMesh.tag_get_handle("psi", 1, MB_TYPE_DOUBLE, th_psi,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
Tag th_rho;
CHKERR postProcMesh.tag_get_handle("RHO", 1, MB_TYPE_DOUBLE, th_rho,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
Tag th_piola1;
Tag principal;
Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
Tag th_R;
Tag th_grad_rho;
if (!commonData.less_post_proc) {
CHKERR postProcMesh.tag_get_handle("Piola_Stress1", 9, MB_TYPE_DOUBLE,
th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
CHKERR postProcMesh.tag_get_handle("Principal_stress", 3, MB_TYPE_DOUBLE,
principal, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
CHKERR postProcMesh.tag_get_handle("Principal_stress_S1", 3, MB_TYPE_DOUBLE,
th_prin_stress_vect1,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("Principal_stress_S2", 3, MB_TYPE_DOUBLE,
th_prin_stress_vect2,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("Principal_stress_S3", 3, MB_TYPE_DOUBLE,
th_prin_stress_vect3,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("R", 1, MB_TYPE_DOUBLE, th_R,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("grad_rho", 3, MB_TYPE_DOUBLE,
th_grad_rho,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
}
const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
MatrixDouble &matS = commonData.data.matS;
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
auto psi = getFTensor0FromVec(commonData.data.vecPsi);
auto R = getFTensor0FromVec(commonData.data.vecR);
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
MatrixDouble3by3 stress(3, 3);
VectorDouble3 eigen_values(3);
MatrixDouble3by3 eigen_vectors(3, 3);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double val = psi;
CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &val);
val = rho;
CHKERR postProcMesh.tag_set_data(th_rho, &mapGaussPts[gg], 1, &val);
// S is symmetric tensor, need to map it to non-symmetric tensor to save it
// on moab tag to be viable in paraview.
for (int ii = 0; ii != 3; ii++) {
for (int jj = 0; jj != 3; jj++) {
stress(ii, jj) = S(ii, jj);
}
}
CHKERR postProcMesh.tag_set_data(th_piola2, &mapGaussPts[gg], 1,
&stress(0, 0));
if (!commonData.less_post_proc) {
for (int ii = 0; ii != 3; ii++) {
for (int jj = 0; jj != 3; jj++) {
stress(ii, jj) = P(ii, jj);
}
}
CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
&stress(0, 0));
val = R;
CHKERR postProcMesh.tag_set_data(th_R, &mapGaussPts[gg], 1, &val);
const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
CHKERR postProcMesh.tag_set_data(th_grad_rho, &mapGaussPts[gg], 1, t2);
// Add calculation of eigen values, i.e. prinple stresses.
noalias(eigen_vectors) = stress;
VectorDouble3 prin_stress_vect1(3);
VectorDouble3 prin_stress_vect2(3);
VectorDouble3 prin_stress_vect3(3);
// LAPACK - eigenvalues and vectors. Applied twice for initial creates
// memory space
int n = 3, lda = 3, info, lwork = -1;
double wkopt;
info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
&(eigen_values.data()[0]), &wkopt, lwork);
if (info != 0)
SETERRQ1(PETSC_COMM_SELF, 1,
"something is wrong with lapack_dsyev info = %d", info);
lwork = (int)wkopt;
double work[lwork];
info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
&(eigen_values.data()[0]), work, lwork);
if (info != 0)
SETERRQ1(PETSC_COMM_SELF, 1,
"something is wrong with lapack_dsyev info = %d", info);
for (int ii = 0; ii < 3; ii++) {
prin_stress_vect1[ii] = eigen_vectors(0, ii);
prin_stress_vect2[ii] = eigen_vectors(1, ii);
prin_stress_vect3[ii] = eigen_vectors(2, ii);
}
CHKERR postProcMesh.tag_set_data(principal, &mapGaussPts[gg], 1,
&eigen_values[0]);
CHKERR postProcMesh.tag_set_data(th_prin_stress_vect1, &mapGaussPts[gg],
1, &*prin_stress_vect1.begin());
CHKERR postProcMesh.tag_set_data(th_prin_stress_vect2, &mapGaussPts[gg],
1, &*prin_stress_vect2.begin());
CHKERR postProcMesh.tag_set_data(th_prin_stress_vect3, &mapGaussPts[gg],
1, &*prin_stress_vect3.begin());
}
++S;
++P;
++psi;
++R;
++rho;
++grad_rho;
}
}
OpCalculateStressTangent::OpCalculateStressTangent(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {
// Set that this operator is executed only for vertices on element
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
OpCalculateStressTangent::doWork(int side, EntityType type,
if (type != MBVERTEX)
auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
commonData.data.matInvF.resize(9, nb_gauss_pts, false);
auto invF = getFTensor2FromMat<3, 3>(commonData.data.matInvF);
commonData.data.vecDetF.resize(nb_gauss_pts, false);
auto det_f = getFTensor0FromVec(commonData.data.vecDetF);
MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
dP_dF.resize(81, nb_gauss_pts, false);
MatrixDouble &matF = commonData.data.matGradientOfDeformation;
matF.resize(9, nb_gauss_pts, false);
commonData.data.matP.resize(9, nb_gauss_pts, false);
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
commonData.data.vecPsi.resize(nb_gauss_pts, false);
auto psi = getFTensor0FromVec(commonData.data.vecPsi);
const double lambda = commonData.lambda;
const double mu = commonData.mu;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get deformation gradient
F(i, I) = grad(i, I);
F(0, 0) += 1;
F(1, 1) += 1;
F(2, 2) += 1;
CHKERR invertTensor3by3(F, det_f, invF);
const double log_det = log(det_f);
psi = F(i, J) * F(i, J) - 3 - 2 * log_det;
psi *= 0.5 * mu;
psi += 0.5 * lambda * log_det * log_det;
const double var = lambda * log_det - mu;
P(i, J) = mu * F(i, J) + var * invF(J, i);
const double coef = mu - lambda * log_det;
// TODO: This can be improved by utilising the symmetries of the D2 tensor
D2(i, J, k, L) =
invF(J, i) * (invF(L, k) * lambda) + invF(L, i) * (invF(J, k) * coef);
D2(0, 0, 0, 0) += mu;
D2(0, 1, 0, 1) += mu;
D2(0, 2, 0, 2) += mu;
D2(1, 0, 1, 0) += mu;
D2(1, 1, 1, 1) += mu;
D2(1, 2, 1, 2) += mu;
D2(2, 0, 2, 0) += mu;
D2(2, 1, 2, 1) += mu;
D2(2, 2, 2, 2) += mu;
++det_f;
++invF;
++grad;
++F;
++D2;
++P;
++psi;
}
}
OpAssmbleStressRhs::OpAssmbleStressRhs(Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {}
OpAssmbleStressRhs::doWork(int side, EntityType type,
const int nb_dofs = data.getIndices().size();
if (!nb_dofs)
nF.resize(nb_dofs, false);
nF.clear();
const int nb_gauss_pts = data.getN().size1();
const int nb_base_functions = data.getN().size2();
auto diff_base_functions = data.getFTensor1DiffN<3>();
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
// MatrixDouble &matS = commonData.data.matS;
// MatrixDouble &matF = commonData.data.matGradientOfDeformation;
// FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
// FTensor::Tensor2_symmetric<double*,3> S(SYMMETRIC_TENSOR2_MAT_PTR(matS),6);
const double n = commonData.n;
const double rho_ref = commonData.rHo_ref;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
}
double a = w * pow(rho / rho_ref, n);
int bb = 0;
for (; bb != nb_dofs / 3; bb++) {
double *ptr = &nF[3 * bb];
FTensor::Tensor1<double *, 3> t1(ptr, &ptr[1], &ptr[2]);
t1(i) += a * P(i, I) * diff_base_functions(I);
// t1(i) += 0.5*a*diff_base_functions(J)*(F(i,I)*S(I,J));
// t1(i) += 0.5*a*diff_base_functions(I)*(F(i,J)*S(I,J));
++diff_base_functions;
}
// Could be that we base more base functions than needed to approx. disp.
// field.
for (; bb != nb_base_functions; bb++) {
++diff_base_functions;
}
++P;
// ++S;
// ++F;
++rho;
}
// Assemble internal force vector for time solver (TS)
getFEMethod()->ts_F, /// Get the right hand side vector for time solver
nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
}
OpAssmbleRhoRhs::OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
"RHO", UserDataOperator::OPROW),
commonData(common_data) {}
OpAssmbleRhoRhs::doWork(int side, EntityType type,
const int nb_dofs = data.getIndices().size();
if (!nb_dofs)
nF.resize(nb_dofs, false);
nF.clear();
const int nb_gauss_pts = data.getN().size1();
const int nb_base_functions = data.getN().size2();
auto base_functions = data.getFTensor0N();
auto diff_base_functions = data.getFTensor1DiffN<3>();
if (!commonData.data.rhoPtr) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "rhoPtsr is null");
}
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
if (!commonData.data.gradRhoPtr) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "gradRhoPtr is null");
}
auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
auto R = getFTensor0FromVec(commonData.data.vecR);
if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
commonData.data.vecRhoDt.clear();
// SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Time derivative of
// density at integration not calculated");
}
auto drho_dt = getFTensor0FromVec(commonData.data.vecRhoDt);
const double R0 = commonData.R0;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
}
// cerr << R << endl;
int bb = 0;
for (; bb != nb_dofs; bb++) {
nF[bb] += w * base_functions * drho_dt; // Density rate term
nF[bb] -= w * base_functions * R; // Source term
nF[bb] -= w * R0 * diff_base_functions(I) * grad_rho(I); // Gradient term
++base_functions;
++diff_base_functions;
}
// Could be that we base more base functions than needed to approx. density
// field.
for (; bb != nb_base_functions; bb++) {
++base_functions;
++diff_base_functions;
}
// Increment quantities for integration pts.
++rho;
++grad_rho;
++drho_dt;
++R;
}
// Assemble internal force vector for time solver (TS)
getFEMethod()->ts_F, /// Get the right hand side vector for time solver
nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
}
OpAssmbleRhoLhs_dRho::OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
"RHO", "RHO", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = true;
}
OpAssmbleRhoLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
// Set size can clear local tangent matrix
locK_rho_rho.resize(row_nb_dofs, col_nb_dofs, false);
locK_rho_rho.clear();
const int row_nb_gauss_pts = row_data.getN().size1();
if (!row_nb_gauss_pts)
const int row_nb_base_functions = row_data.getN().size2();
const double ts_a = getFEMethod()->ts_a;
const double R0 = commonData.R0;
const double c = commonData.c;
const double n = commonData.n;
const double m = commonData.m;
const double rho_ref = commonData.rHo_ref;
const double psi_ref = commonData.pSi_ref;
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
auto psi = getFTensor0FromVec(commonData.data.vecPsi);
auto row_base_functions = row_data.getFTensor0N();
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
}
const double kp = commonData.getCFromDensity(rho);
const double kp_diff = commonData.getCFromDensityDiff(rho);
const double a = c * kp * ((n - m) / rho) * pow(rho / rho_ref, n - m) * psi;
const double a_diff =
c * kp_diff * (pow(rho / rho_ref, n - m) * psi - psi_ref);
int row_bb = 0;
for (; row_bb != row_nb_dofs; row_bb++) {
auto col_base_functions = col_data.getFTensor0N(gg, 0);
auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
locK_rho_rho(row_bb, col_bb) +=
ts_a * w * row_base_functions * col_base_functions;
locK_rho_rho(row_bb, col_bb) -=
w * R0 * row_diff_base_functions(I) * col_diff_base_functions(I);
locK_rho_rho(row_bb, col_bb) -=
a * w * row_base_functions * col_base_functions;
locK_rho_rho(row_bb, col_bb) -=
a_diff * w * row_base_functions * col_base_functions;
++col_base_functions;
++col_diff_base_functions;
}
++row_base_functions;
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_base_functions;
++row_diff_base_functions;
}
++psi;
++rho;
}
// cerr << locK_rho_rho << endl;
CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
&*row_data.getIndices().begin(), col_nb_dofs,
&*col_data.getIndices().begin(), &locK_rho_rho(0, 0),
ADD_VALUES);
// is symmetric
if (row_side != col_side || row_type != col_type) {
transLocK_rho_rho.resize(col_nb_dofs, row_nb_dofs, false);
noalias(transLocK_rho_rho) = trans(locK_rho_rho);
CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs,
&*col_data.getIndices().begin(), row_nb_dofs,
&*row_data.getIndices().begin(),
&transLocK_rho_rho(0, 0), ADD_VALUES);
}
}
OpAssmbleRhoLhs_dF::OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
"RHO", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = false;
}
OpAssmbleRhoLhs_dF::doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
// Set size can clear local tangent matrix
locK_rho_F.resize(row_nb_dofs, col_nb_dofs, false);
locK_rho_F.clear();
const int row_nb_gauss_pts = row_data.getN().size1();
const int col_nb_base_functions = col_data.getN().size2();
const double c = commonData.c;
const double n = commonData.n;
const double m = commonData.m;
const double rho_ref = commonData.rHo_ref;
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
}
const double kp = commonData.getCFromDensity(rho);
const double a = w * c * kp * pow(rho / rho_ref, n - m);
int col_bb = 0;
for (; col_bb != col_nb_dofs / 3; col_bb++) {
t1(i) = a * P(i, I) * col_diff_base_functions(I);
auto row_base_functions = row_data.getFTensor0N(gg, 0);
for (int row_bb = 0; row_bb != row_nb_dofs; row_bb++) {
FTensor::Tensor1<double *, 3> k(&locK_rho_F(row_bb, 3 * col_bb + 0),
&locK_rho_F(row_bb, 3 * col_bb + 1),
&locK_rho_F(row_bb, 3 * col_bb + 2));
k(i) -= row_base_functions * t1(i);
++row_base_functions;
}
++col_diff_base_functions;
}
for (; col_bb != col_nb_base_functions; col_bb++) {
++col_diff_base_functions;
}
++P;
++rho;
}
CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
&*row_data.getIndices().begin(), col_nb_dofs,
&*col_data.getIndices().begin(), &locK_rho_F(0, 0),
ADD_VALUES);
}
template <bool ADOLC>
OpAssmbleStressLhs_dF<ADOLC>::OpAssmbleStressLhs_dF(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = true;
}
template <bool ADOLC>
MoFEMErrorCode OpAssmbleStressLhs_dF<ADOLC>::doWork(
int row_side, int col_side, EntityType row_type, EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
const bool diagonal_block = (row_type == col_type) && (row_side == col_side);
// Set size can clear local tangent matrix
locK_P_F.resize(row_nb_dofs, col_nb_dofs, false);
locK_P_F.clear();
const int row_nb_gauss_pts = row_data.getN().size1();
const int row_nb_base_functions = row_data.getN().size2();
MatrixDouble &matS = commonData.data.matS;
MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
// MatrixDouble &dS_dC = commonData.data.matMaterialTangent;
// MatrixDouble &matF = commonData.data.matGradientOfDeformation;
double t0;
const double n = commonData.n;
const double rho_ref = commonData.rHo_ref;
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
// FTensor::Tensor4_ddg<double*,3,3> D1(SYMMETRIC_TENSOR4_MAT_PTR(dS_dC),6);
// FTensor::Tensor2<double*,3,3> F(TENSOR2_MAT_PTR(matF));
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
}
const double a = w * pow(rho / rho_ref, n);
int row_bb = 0;
for (; row_bb != row_nb_dofs / 3; row_bb++) {
auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
int col_bb = 0;
for (; col_bb != final_bb; col_bb++) {
&locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
&locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
&locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
&locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
&locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
&locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
&locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
&locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
&locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
diffDiff(J, L) =
a * row_diff_base_functions(J) * col_diff_base_functions(L);
// TODO: FIXME: Tensor D2 has major symmetry, we do not have such tensor
// implemented at the moment. Using tensor with major symmetry will
// improve code efficiency.
t_assemble(i, k) += diffDiff(J, L) * D2(i, J, k, L);
if (ADOLC) {
// Stress stiffening part (only with adolc since it is using dS / dC
// derrivative. Check OpCalculateStressTangentWithAdolc operator)
t0 = diffDiff(J, L) * S(J, L);
t_assemble(0, 0) += t0;
t_assemble(1, 1) += t0;
t_assemble(2, 2) += t0;
}
// Next base function for column
++col_diff_base_functions;
}
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_diff_base_functions;
}
// ++D1;
++D2;
++S;
// ++F;
++rho;
}
// Copy of symmetric part for assemble
if (diagonal_block) {
for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
int col_bb = 0;
for (; col_bb != row_bb + 1; col_bb++) {
&locK_P_F(3 * row_bb + 0, 3 * col_bb + 0),
&locK_P_F(3 * row_bb + 0, 3 * col_bb + 1),
&locK_P_F(3 * row_bb + 0, 3 * col_bb + 2),
&locK_P_F(3 * row_bb + 1, 3 * col_bb + 0),
&locK_P_F(3 * row_bb + 1, 3 * col_bb + 1),
&locK_P_F(3 * row_bb + 1, 3 * col_bb + 2),
&locK_P_F(3 * row_bb + 2, 3 * col_bb + 0),
&locK_P_F(3 * row_bb + 2, 3 * col_bb + 1),
&locK_P_F(3 * row_bb + 2, 3 * col_bb + 2));
&locK_P_F(3 * col_bb + 0, 3 * row_bb + 0),
&locK_P_F(3 * col_bb + 0, 3 * row_bb + 1),
&locK_P_F(3 * col_bb + 0, 3 * row_bb + 2),
&locK_P_F(3 * col_bb + 1, 3 * row_bb + 0),
&locK_P_F(3 * col_bb + 1, 3 * row_bb + 1),
&locK_P_F(3 * col_bb + 1, 3 * row_bb + 2),
&locK_P_F(3 * col_bb + 2, 3 * row_bb + 0),
&locK_P_F(3 * col_bb + 2, 3 * row_bb + 1),
&locK_P_F(3 * col_bb + 2, 3 * row_bb + 2));
t_off_side(i, k) = t_assemble(k, i);
}
}
}
const int *row_ind = &*row_data.getIndices().begin();
const int *col_ind = &*col_data.getIndices().begin();
CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs, row_ind, col_nb_dofs,
col_ind, &locK_P_F(0, 0), ADD_VALUES);
if (row_type != col_type || row_side != col_side) {
transLocK_P_F.resize(col_nb_dofs, row_nb_dofs, false);
noalias(transLocK_P_F) = trans(locK_P_F);
CHKERR MatSetValues(getFEMethod()->ts_B, col_nb_dofs, col_ind, row_nb_dofs,
row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
}
}
OpAssmbleStressLhs_dRho::OpAssmbleStressLhs_dRho(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", "RHO", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = false; // This is off-diagonal (upper-left), so no symmetric
}
OpAssmbleStressLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
// Set size can clear local tangent matrix
locK_P_RHO.resize(row_nb_dofs, col_nb_dofs, false);
locK_P_RHO.clear();
const double n = commonData.n;
const double rho_ref = commonData.rHo_ref;
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
const int row_nb_gauss_pts = row_data.getN().size1();
const int row_nb_base_functions = row_data.getN().size2();
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
}
double a = w * (n / rho) * pow(rho / rho_ref, n);
int row_bb = 0;
for (; row_bb != row_nb_dofs / 3; row_bb++) {
t1(i) = a * row_diff_base_functions(I) * P(i, I);
auto base_functions = col_data.getFTensor0N(gg, 0);
for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
// if(base_functions!=col_data.getN(gg)[col_bb]) {
// cerr << "Error" << endl;
// }
FTensor::Tensor1<double *, 3> k(&locK_P_RHO(3 * row_bb + 0, col_bb),
&locK_P_RHO(3 * row_bb + 1, col_bb),
&locK_P_RHO(3 * row_bb + 2, col_bb));
k(i) += t1(i) * base_functions;
++base_functions;
}
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_diff_base_functions;
}
++P;
++rho;
}
CHKERR MatSetValues(getFEMethod()->ts_B, row_nb_dofs,
&*row_data.getIndices().begin(), col_nb_dofs,
&*col_data.getIndices().begin(), &locK_P_RHO(0, 0),
ADD_VALUES);
}
Remodeling::CommonData &common_data)
: mField(m_field), commonData(common_data), postProc(m_field),
postProcElastic(m_field),
// densityMaps(m_field),
iNit(false) {}
}
}
PetscBool mass_postproc = PETSC_FALSE;
PetscBool equilibrium_flg = PETSC_FALSE;
PetscBool analysis_mesh_flg = PETSC_FALSE;
rate = 1;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling post-process",
"none");
CHKERR PetscOptionsBool(
"-mass_postproc",
"is used for testing, calculates mass and energy at each time step", "",
mass_postproc, &mass_postproc, PETSC_NULL);
CHKERR PetscOptionsBool("-analysis_mesh",
"saves analysis mesh at each time step", "",
analysis_mesh_flg, &analysis_mesh_flg, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-equilibrium_stop_rate",
"is used to stop calculations when equilibium state is achieved", "",
rate, &rate, &equilibrium_flg);
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
if (!iNit) {
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
"Bone remodeling post-process", "none");
pRT = 1;
CHKERR PetscOptionsInt("-my_output_prt",
"frequncy how often results are dumped on hard disk",
"", 1, &pRT, NULL);
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
CHKERR postProc.generateReferenceElementMesh();
CHKERR postProc.addFieldValuesPostProc("DISPLACEMENTS");
CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
if (!commonData.less_post_proc) {
CHKERR postProc.addFieldValuesGradientPostProc("DISPLACEMENTS");
}
// CHKERR postProc.addFieldValuesPostProc("RHO");
// CHKERR postProc.addFieldValuesGradientPostProc("RHO");
{
boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
&list_of_operators = postProc.getOpPtrVector();
list_of_operators.push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
list_of_operators.push_back(new OpCalculateScalarFieldGradient<3>(
"RHO", commonData.data.gradRhoPtr));
// Get displacement gradient at integration points
list_of_operators.push_back(new OpCalculateVectorFieldGradient<3, 3>(
"DISPLACEMENTS", commonData.data.gradDispPtr));
list_of_operators.push_back(new OpCalculateStress(commonData));
list_of_operators.push_back(new OpPostProcStress(
postProc.postProcMesh, postProc.mapGaussPts, commonData));
}
CHKERR postProcElastic.generateReferenceElementMesh();
CHKERR postProcElastic.addFieldValuesPostProc("DISPLACEMENTS");
CHKERR postProcElastic.addFieldValuesPostProc("MESH_NODE_POSITIONS");
CHKERR postProcElastic.addFieldValuesGradientPostProc("DISPLACEMENTS");
std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
commonData.elasticPtr->setOfBlocks.begin();
for (; sit != commonData.elasticPtr->setOfBlocks.end(); sit++) {
postProcElastic.getOpPtrVector().push_back(new PostProcStress(
postProcElastic.postProcMesh, postProcElastic.mapGaussPts,
"DISPLACEMENTS", sit->second, postProcElastic.commonData));
}
iNit = true;
}
if (mass_postproc) {
Vec mass_vec;
Vec energ_vec;
int mass_vec_ghost[] = {0};
CHKERR VecCreateGhost(PETSC_COMM_WORLD, (!mField.get_comm_rank()) ? 1 : 0,
1, 1, mass_vec_ghost, &mass_vec);
CHKERR VecZeroEntries(mass_vec);
CHKERR VecDuplicate(mass_vec, &energ_vec);
Remodeling::Fe energyCalc(mField);
energyCalc.getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
"RHO", commonData.data.gradRhoPtr));
energyCalc.getOpPtrVector().push_back(
commonData.data.gradDispPtr));
energyCalc.getOpPtrVector().push_back(new OpCalculateStress(commonData));
energyCalc.getOpPtrVector().push_back(
new OpMassAndEnergyCalculation("RHO", commonData, energ_vec, mass_vec));
CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING",
&energyCalc);
CHKERR VecAssemblyBegin(mass_vec);
CHKERR VecAssemblyEnd(mass_vec);
double mass_sum;
CHKERR VecSum(mass_vec, &mass_sum);
CHKERR VecAssemblyBegin(energ_vec);
CHKERR VecAssemblyEnd(energ_vec);
double energ_sum;
CHKERR VecSum(energ_vec, &energ_sum);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Time: %g Mass: %6.5g Elastic energy: %6.5g \n", ts_t,
mass_sum, energ_sum);
CHKERR VecDestroy(&mass_vec);
CHKERR VecDestroy(&energ_vec);
if (equilibrium_flg) {
double equilibrium_rate =
fabs(energ_sum - commonData.cUrrent_psi) / energ_sum;
double equilibrium_mass_rate =
fabs(mass_sum - commonData.cUrrent_mass) / mass_sum;
if (equilibrium_rate < rate) {
CHKERR PetscPrintf(
PETSC_COMM_WORLD,
"Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
equilibrium_rate * 100);
}
if (equilibrium_mass_rate < rate) {
CHKERR PetscPrintf(
PETSC_COMM_WORLD,
"Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
equilibrium_mass_rate * 100);
}
commonData.cUrrent_psi = energ_sum;
commonData.cUrrent_mass = mass_sum;
}
}
int step;
#if PETSC_VERSION_LE(3, 8, 0)
CHKERR TSGetTimeStepNumber(ts, &step);
#else
CHKERR TSGetStepNumber(ts, &step);
#endif
if ((step) % pRT == 0) {
ostringstream sss;
sss << "out_" << step << ".h5m";
CHKERR DMoFEMLoopFiniteElements(commonData.dm, "FE_REMODELLING", &postProc);
CHKERR postProc.postProcMesh.write_file(sss.str().c_str(), "MOAB",
"PARALLEL=WRITE_PART");
if (analysis_mesh_flg) {
ostringstream ttt;
ttt << "analysis_mesh_" << step << ".h5m";
CHKERR mField.get_moab().write_file(ttt.str().c_str(), "MOAB",
"PARALLEL=WRITE_PART");
}
if (commonData.nOremodellingBlock) {
CHKERR DMoFEMLoopFiniteElements(commonData.dm, "ELASTIC",
&postProcElastic);
ostringstream ss;
ss << "out_elastic_" << step << ".h5m";
CHKERR postProcElastic.postProcMesh.write_file(ss.str().c_str(), "MOAB",
"PARALLEL=WRITE_PART");
}
}
}
MoFEMErrorCode Remodeling::getParameters() {
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
commonData.oRder = 2;
CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
&commonData.oRder, PETSC_NULL);
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
string name = it->getName();
if (name.compare(0, 14, "NO_REMODELLING") == 0) {
commonData.nOremodellingBlock = true;
EntityHandle meshset = it->getMeshset();
CHKERR this->mField.get_moab().get_entities_by_type(
meshset, MBTET, commonData.tEts_block, true);
commonData.tEts_all =
subtract(commonData.tEts_all, commonData.tEts_block);
}
}
}
MoFEMErrorCode Remodeling::addFields() {
// Seed all mesh entities to MoFEM database, those entities can be potentially
// used as finite elements or as entities which carry some approximation
// field.
commonData.bitLevel.set(0);
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 3, commonData.bitLevel);
int order = commonData.oRder;
// Add displacement field
CHKERR mField.add_field("DISPLACEMENTS", H1,
/*AINSWORTH_LOBATTO_BASE*/ AINSWORTH_LEGENDRE_BASE, 3,
MB_TAG_SPARSE, MF_ZERO);
// Add field representing ho-geometry
CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
MB_TAG_SPARSE, MF_ZERO);
// Check if density is available, if not add density field.
bool add_rho_field = false;
if (!mField.check_field("RHO")) {
CHKERR mField.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1,
MB_TAG_SPARSE, MF_ZERO);
// FIXME
CHKERR mField.add_ents_to_field_by_type(commonData.tEts_all, MBTET, "RHO");
// CHKERR mField.add_ents_to_field_by_type(0,MBTET,"RHO");
CHKERR mField.synchronise_field_entities("RHO");
add_rho_field = true;
CHKERR mField.set_field_order(0, MBVERTEX, "RHO", 1);
CHKERR mField.set_field_order(0, MBEDGE, "RHO", order - 1);
CHKERR mField.set_field_order(0, MBTRI, "RHO", order - 1);
CHKERR mField.set_field_order(0, MBTET, "RHO", order - 1);
}
// Add entities to field
CHKERR mField.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENTS");
CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
// Set approximation order to entities
CHKERR mField.set_field_order(0, MBVERTEX, "DISPLACEMENTS", 1);
CHKERR mField.set_field_order(0, MBEDGE, "DISPLACEMENTS", order);
CHKERR mField.set_field_order(0, MBTRI, "DISPLACEMENTS", order);
CHKERR mField.set_field_order(0, MBTET, "DISPLACEMENTS", order);
// Assumes that geometry is approximated using 2nd order polynomials.
// Apply 2nd order only on skin
{
// Skinner skin(&mField.get_moab());
// Range faces,tets;
// CHKERR mField.get_moab().get_entities_by_type(0,MBTET,tets);
// CHKERR skin.find_skin(0,tets,false,faces);
// Range edges;
// CHKERR mField.get_moab().get_adjacencies(
// faces,1,false,edges,moab::Interface::UNION
// );
CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
}
CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
// Build fields
CHKERR mField.build_fields();
// If order was not set from CT scans set homogenous order equal to
// reference bone density
if (add_rho_field) {
CHKERR mField.getInterface<FieldBlas>()->setField(commonData.rHo_ref,
MBVERTEX, "RHO");
// const DofEntity_multiIndex *dofs_ptr;
// CHKERR mField.get_dofs(&dofs_ptr);
// for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(mField,"RHO",dit)) {
// cerr << (*dit)->getFieldData() << endl;
// }
}
// Project geometry given on 10-node tets on ho-geometry
Projection10NodeCoordsOnField ent_method_material(mField,
"MESH_NODE_POSITIONS");
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
}
MoFEMErrorCode Remodeling::addElements() {
CHKERR mField.add_finite_element("FE_REMODELLING", MF_ZERO);
CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING", "RHO");
CHKERR mField.modify_finite_element_add_field_data("FE_REMODELLING",
"MESH_NODE_POSITIONS");
CHKERR mField.modify_finite_element_add_field_row("FE_REMODELLING", "RHO");
CHKERR mField.modify_finite_element_add_field_col("FE_REMODELLING", "RHO");
CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_all, MBTET,
"FE_REMODELLING");
CHKERR mField.add_finite_element("ELASTIC");
CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
"MESH_NODE_POSITIONS");
if (commonData.nOremodellingBlock) {
CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_block,
MBTET, "ELASTIC");
}
commonData.elasticMaterialsPtr =
boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
new NonlinearElasticElement(mField, 10));
CHKERR commonData.elasticMaterialsPtr->setBlocks(
commonData.elasticPtr->setOfBlocks);
CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
CHKERR commonData.elasticPtr->setOperators(
"DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
CHKERR mField.build_finite_elements();
CHKERR mField.build_adjacencies(commonData.bitLevel);
// Allocate memory for density and gradient of displacements at integration
// points
commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
commonData.data.gradRhoPtr =
boost::shared_ptr<MatrixDouble>(new MatrixDouble());
commonData.data.gradDispPtr =
boost::shared_ptr<MatrixDouble>(new MatrixDouble());
// Add operators Rhs
boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
&list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
// Get density at integration points
list_of_operators_rhs.push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
list_of_operators_rhs.push_back(
new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
// Get displacement gradient at integration points
list_of_operators_rhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
"DISPLACEMENTS", commonData.data.gradDispPtr));
list_of_operators_rhs.push_back(new OpCalculateStress(commonData));
list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
// Add operators Lhs
boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
&list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
// Get density at integration points
list_of_operators_lhs.push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
list_of_operators_rhs.push_back(
new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
// Get displacement gradient at integration points
list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
"DISPLACEMENTS", commonData.data.gradDispPtr));
if (commonData.with_adol_c) {
list_of_operators_lhs.push_back(
new OpCalculateStressTangentWithAdolc(commonData));
list_of_operators_lhs.push_back(
new OpAssmbleStressLhs_dF<true>(commonData));
} else {
list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
list_of_operators_lhs.push_back(
new OpAssmbleStressLhs_dF<false>(commonData));
}
list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
}
MoFEMErrorCode Remodeling::addMomentumFluxes() {
// Add Neumann forces elements
CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
CHKERR MetaNodalForces::addElement(mField, "DISPLACEMENTS");
CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
// forces and pressures on surface
CHKERR MetaNeummanForces::setMomentumFluxOperators(
mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
// noadl forces
CHKERR MetaNodalForces::setOperators(mField, commonData.nodalForces,
PETSC_NULL, "DISPLACEMENTS");
// edge forces
CHKERR MetaEdgeForces::setOperators(mField, commonData.edgeForces, PETSC_NULL,
"DISPLACEMENTS");
for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
commonData.neumannForces.begin();
mit != commonData.neumannForces.end(); mit++) {
mit->second->methodsOp.push_back(
new TimeForceScale("-my_load_history", false));
}
for (boost::ptr_map<string, NodalForce>::iterator mit =
commonData.nodalForces.begin();
mit != commonData.nodalForces.end(); mit++) {
mit->second->methodsOp.push_back(
new TimeForceScale("-my_load_history", false));
}
for (boost::ptr_map<string, EdgeForce>::iterator mit =
commonData.edgeForces.begin();
mit != commonData.edgeForces.end(); mit++) {
mit->second->methodsOp.push_back(
new TimeForceScale("-my_load_history", false));
}
}
MoFEMErrorCode Remodeling::buildDM() {
commonData.dm_name = "DM_REMODELING";
// Register DM problem
CHKERR DMRegister_MoFEM(commonData.dm_name);
CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
CHKERR DMSetType(commonData.dm, commonData.dm_name);
CHKERR DMMoFEMSetIsPartitioned(commonData.dm, PETSC_TRUE);
// Create DM instance
CHKERR DMMoFEMCreateMoFEM(commonData.dm, &mField, commonData.dm_name,
commonData.bitLevel);
// Configure DM form line command options (DM itself, solvers,
// pre-conditioners, ... )
CHKERR DMSetFromOptions(commonData.dm);
// Add elements to dm (only one here)
CHKERR DMMoFEMAddElement(commonData.dm, "FE_REMODELLING");
CHKERR DMMoFEMAddElement(commonData.dm, "ELASTIC");
if (mField.check_finite_element("FORCE_FE")) {
CHKERR DMMoFEMAddElement(commonData.dm, "FORCE_FE");
}
if (mField.check_finite_element("PRESSURE_FE")) {
CHKERR DMMoFEMAddElement(commonData.dm, "PRESSURE_FE");
}
mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
CHKERR DMSetUp(commonData.dm);
}
MoFEMErrorCode Remodeling::solveDM() {
Vec F = commonData.F;
Vec D = commonData.D;
Mat A = commonData.A;
DM dm = commonData.dm;
TS ts = commonData.ts;
CHKERR TSSetIFunction(ts, F, PETSC_NULL, PETSC_NULL);
CHKERR TSSetIJacobian(ts, A, A, PETSC_NULL, PETSC_NULL);
double ftime = 1;
#if PETSC_VERSION_GE(3, 8, 0)
CHKERR TSSetMaxTime(ts, ftime);
#endif
CHKERR TSSetFromOptions(ts);
CHKERR TSSetDM(ts, dm);
SNES snes;
CHKERR TSGetSNES(ts, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Set up FIELDSPLIT
// Only is user set -pc_type fieldsplit
if (is_pcfs == PETSC_TRUE) {
IS is_disp, is_rho;
const MoFEM::Problem *problem_ptr;
CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
problem_ptr->getName(), ROW, "DISPLACEMENTS", 0, 3, &is_disp);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
problem_ptr->getName(), ROW, "RHO", 0, 1, &is_rho);
CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
CHKERR ISDestroy(&is_disp);
CHKERR ISDestroy(&is_rho);
}
// Monitor
MonitorPostProc post_proc(mField, commonData);
TsCtx *ts_ctx;
DMMoFEMGetTsCtx(dm, &ts_ctx);
{
ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
CHKERR TSMonitorSet(ts, f_TSMonitorSet, ts_ctx, PETSC_NULL);
}
#if PETSC_VERSION_GE(3, 7, 0)
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
CHKERR TSSolve(ts, D);
CHKERR TSGetTime(ts, &ftime);
PetscInt steps, snesfails, rejects, nonlinits, linits;
#if PETSC_VERSION_LE(3, 8, 0)
CHKERR TSGetTimeStepNumber(ts, &steps);
#else
CHKERR TSGetStepNumber(ts, &steps);
#endif
CHKERR TSGetSNESFailures(ts, &snesfails);
CHKERR TSGetStepRejections(ts, &rejects);
CHKERR TSGetSNESIterations(ts, &nonlinits);
CHKERR TSGetKSPIterations(ts, &linits);
PetscPrintf(PETSC_COMM_WORLD,
"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
"linits %D\n",
steps, rejects, snesfails, ftime, nonlinits, linits);
if (commonData.is_atom_testing) {
if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
}
}
MoFEMErrorCode OpMassAndEnergyCalculation::doWork(
int row_side, EntityType row_type,
if (row_type != MBVERTEX)
const int nb_gauss_pts = row_data.getN().size1();
// commonData.data.vecPsi.resize(nb_gauss_pts,false);
auto rho = getFTensor0FromVec(*commonData.data.rhoPtr);
auto psi = getFTensor0FromVec(commonData.data.vecPsi);
double energy = 0;
double mass = 0;
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double vol = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
vol *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
}
energy += vol * psi;
mass += vol * rho;
++psi;
++rho;
}
CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
}
} // namespace BoneRemodeling
static Index< 'm', 3 > m
static Index< 'n', 3 > n
static Index< 'L', 3 > L
static Index< 'J', 3 > J
static Index< 'I', 3 > I
static Index< 'K', 3 > K
Elastic materials.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
Definition: Remodeling.cpp:57
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
Definition: Remodeling.cpp:54
#define TENSOR2_MAT_PTR(MAT)
Definition: Remodeling.cpp:50
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
Definition: Remodeling.cpp:40
#define TENSOR4_MAT_PTR(MAT)
Definition: Remodeling.cpp:48
@ ROW
Definition: definitions.h:192
@ MF_ZERO
Definition: definitions.h:189
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ BLOCKSET
Definition: definitions.h:217
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:131
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:126
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:132
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
#define MU(E, NU)
Definition: fem_tools.h:33
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:982
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:105
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1001
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:886
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:143
DEPRECATED PetscErrorCode f_TSMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Definition: TsCtx.hpp:327
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:905
double w(const double g, const double t)
constexpr auto VecSetValues
constexpr auto MatSetValues
const double D
diffusivity
double young_modulus
Definition: plastic.cpp:103
double rho
Definition: plastic.cpp:105
double poisson_ratio
Definition: plastic.cpp:104
DataForcesAndSourcesCore::EntData EntData
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: EdgeForce.hpp:109
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:74
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:140
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:104
Deprecated interface functions.
keeps basic data about problem
std::string getName() const
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MonitorPostProc(MoFEM::Interface &m_field, std::map< int, NonlinearElasticElement::BlockData > &set_of_blocks, NonlinearElasticElement::MyVolumeFE &fe_elastic_energy, ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, DataAtIntegrationPts &common_data, BlockData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Force scale operator for reading two columns.