v0.8.17
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;
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;
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);
lambda = LAMBDA(young_modulus, poisson_ratio);
mu = MU(young_modulus, poisson_ratio);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Young's modulus E[Pa]: %4.3g\n",
young_modulus);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Poisson ratio nu[-] %4.3g\n",
poisson_ratio);
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();
}
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;
}
}
OpCalulateStressTangentWithAdolc::OpCalulateStressTangentWithAdolc(
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 OpCalulateStressTangentWithAdolc::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;
}
}
OpPostProcStress::OpPostProcStress(moab::Interface &post_proc_mesh,
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;
}
}
OpCalulateStressTangent::OpCalulateStressTangent(
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;
}
OpCalulateStressTangent::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)
CHKERR VecSetValues(
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)
CHKERR VecSetValues(
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 OpCalulateStressTangentWithAdolc 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();
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();
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();
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 OpCalulateStressTangentWithAdolc(commonData));
list_of_operators_lhs.push_back(
new OpAssmbleStressLhs_dF<true>(commonData));
} else {
list_of_operators_lhs.push_back(new OpCalulateStressTangent(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 MetaNodalForces::addElement(mField, "DISPLACEMENTS");
CHKERR MetaEdgeForces::addElement(mField, "DISPLACEMENTS");
// forces and pressures on surface
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