v0.13.1
Remodeling.cpp

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

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 {
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;
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", "",
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));
}
Remodeling::CommonData &common_data)
"RHO", UserDataOperator::OPROW),
commonData(common_data) {}
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) {
}
int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
double *array;
CHKERR VecGetArray(getFEMethod()->ts_u_t, &array);
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);
}
"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;
}
if (type != MBVERTEX)
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();
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Density at integration points is not calculated");
}
commonData.data.vecPsi.resize(nb_gauss_pts, false);
commonData.data.vecR.resize(nb_gauss_pts, false);
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);
matS.resize(nb_gauss_pts, 6, false);
commonData.data.matP.resize(9, nb_gauss_pts, false);
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
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;
}
}
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;
}
if (type != MBVERTEX)
auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
vecC.resize(6, false);
dS_dC.resize(6, 6 * nb_gauss_pts, false);
dS_dC.clear();
matS.resize(nb_gauss_pts, 6, false);
dP_dF.resize(81, nb_gauss_pts, false);
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);
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;
}
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;
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();
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
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));
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;
}
}
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;
}
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);
dP_dF.resize(81, nb_gauss_pts, false);
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);
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;
}
}
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {}
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);
// 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);
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);
}
"RHO", UserDataOperator::OPROW),
commonData(common_data) {}
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>();
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "rhoPtsr is null");
}
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "gradRhoPtr is null");
}
auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
// SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Time derivative of
// density at integration not calculated");
}
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);
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);
}
"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 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);
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);
}
}
"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 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);
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>
Remodeling::CommonData &common_data)
"DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = true;
}
template <bool ADOLC>
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 &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>();
// 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);
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);
}
}
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);
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);
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");
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(
list_of_operators.push_back(new OpCalculateScalarFieldGradient<3>(
// 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);
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", energyCalc, true, false, false,
true);
energyCalc.getOpPtrVector().push_back(
energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
energyCalc.getOpPtrVector().push_back(
energyCalc.getOpPtrVector().push_back(new OpCalculateStress(commonData));
energyCalc.getOpPtrVector().push_back(
new OpMassAndEnergyCalculation("RHO", commonData, energ_vec, mass_vec));
&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;
}
}
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");
}
&postProcElastic);
ostringstream ss;
ss << "out_elastic_" << step << ".h5m";
CHKERR postProcElastic.postProcMesh.write_file(ss.str().c_str(), "MOAB",
"PARALLEL=WRITE_PART");
}
}
}
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Bone remodeling", "none");
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) {
EntityHandle meshset = it->getMeshset();
CHKERR this->mField.get_moab().get_entities_by_type(
meshset, MBTET, commonData.tEts_block, true);
}
}
}
// Seed all mesh entities to MoFEM database, those entities can be potentially
// used as finite elements or as entities which carry some approximation
// field.
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
// 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")) {
MB_TAG_SPARSE, MF_ZERO);
// FIXME
// CHKERR mField.add_ents_to_field_by_type(0,MBTET,"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
// If order was not set from CT scans set homogenous order equal to
// reference bone density
if (add_rho_field) {
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
"MESH_NODE_POSITIONS");
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
}
"DISPLACEMENTS");
"DISPLACEMENTS");
"DISPLACEMENTS");
"MESH_NODE_POSITIONS");
"FE_REMODELLING");
CHKERR mField.modify_finite_element_add_field_row("ELASTIC", "DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_col("ELASTIC", "DISPLACEMENTS");
"DISPLACEMENTS");
"MESH_NODE_POSITIONS");
MBTET, "ELASTIC");
}
boost::shared_ptr<ElasticMaterials>(new ElasticMaterials(mField));
commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
commonData.elasticPtr->setOfBlocks);
CHKERR commonData.elasticPtr->addElement("ELASTIC", "DISPLACEMENTS");
CHKERR commonData.elasticPtr->setOperators(
"DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
// Allocate memory for density and gradient of displacements at integration
// points
commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
boost::shared_ptr<MatrixDouble>(new MatrixDouble());
boost::shared_ptr<MatrixDouble>(new MatrixDouble());
// Add operators Rhs
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feRhs), true, false,
false, false);
auto &list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
// Get density at integration points
list_of_operators_rhs.push_back(
list_of_operators_rhs.push_back(
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
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *(commonData.feLhs), true, false,
false, false);
auto &list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
// Get density at integration points
list_of_operators_lhs.push_back(
list_of_operators_rhs.push_back(
// Get displacement gradient at integration points
list_of_operators_lhs.push_back(new OpCalculateVectorFieldGradient<3, 3>(
"DISPLACEMENTS", commonData.data.gradDispPtr));
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));
}
// Add Neumann forces elements
CHKERR MetaNeummanForces::addNeumannBCElements(mField, "DISPLACEMENTS");
// forces and pressures on surface
CHKERR MetaNeummanForces::setMomentumFluxOperators(
mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
// noadl forces
PETSC_NULL, "DISPLACEMENTS");
// edge forces
"DISPLACEMENTS");
for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
mit != commonData.neumannForces.end(); mit++) {
mit->second->methodsOp.push_back(
new TimeForceScale("-my_load_history", false));
}
for (boost::ptr_map<string, NodalForce>::iterator mit =
mit != commonData.nodalForces.end(); mit++) {
mit->second->methodsOp.push_back(
new TimeForceScale("-my_load_history", false));
}
for (boost::ptr_map<string, EdgeForce>::iterator mit =
mit != commonData.edgeForces.end(); mit++) {
mit->second->methodsOp.push_back(
new TimeForceScale("-my_load_history", false));
}
}
commonData.dm_name = "DM_REMODELING";
// Register DM problem
CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
// Create DM instance
// 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");
if (mField.check_finite_element("FORCE_FE")) {
}
if (mField.check_finite_element("PRESSURE_FE")) {
}
mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
CHKERR DMSetUp(commonData.dm);
}
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 ISSort(is_disp);
CHKERR ISSort(is_rho);
CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
CHKERR ISDestroy(&is_disp);
CHKERR ISDestroy(&is_rho);
}
// Monitor
TsCtx *ts_ctx;
DMMoFEMGetTsCtx(dm, &ts_ctx);
{
ts_ctx->getPostProcessMonitor().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);
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "atom test diverged");
}
}
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);
double energy = 0;
double mass = 0;
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double vol = getVolume() * getGaussPts()(3, gg);
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< '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
constexpr double a
EntitiesFieldData::EntData EntData
@ ROW
Definition: definitions.h:136
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:52
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:47
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define MU(E, NU)
Definition: fem_tools.h:33
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
constexpr double lambda
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1082
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:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:462
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:385
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1101
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:545
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#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.
FTensor::Index< 'i', SPACE_DIM > i
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< '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
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:103
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
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:1218
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DEPRECATED PetscErrorCode f_TSMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Definition: TsCtx.hpp:312
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
double w(const double g, const double t)
OpPlasticTools::CommonData CommonData
const double D
diffusivity
STL namespace.
double A
double young_modulus
Definition: plastic.cpp:105
double rho
Definition: plastic.cpp:107
constexpr double mu
FTensor::Index< 'm', 3 > m
MoFEMErrorCode postProcess()
MoFEMErrorCode operator()()
MoFEMErrorCode preProcess()
MonitorPostProc(MoFEM::Interface &m_field, Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:961
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Remodeling.cpp:969
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:889
OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:883
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleStressLhs_dF(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssmbleStressLhs_dRho(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:820
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:814
OpCalculateStress(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:280
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:294
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:732
OpCalculateStressTangent(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:717
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:396
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:382
OpGetRhoTimeDirevative(Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:232
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:239
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: Remodeling.cpp:567
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
Definition: Remodeling.cpp:550
VectorDouble vecPsi
Elastic energy density.
Definition: Remodeling.hpp:189
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
Definition: Remodeling.hpp:185
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
Definition: Remodeling.hpp:187
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
Definition: Remodeling.hpp:188
MatrixDouble matInvF
inverse of deformation gradient
Definition: Remodeling.hpp:199
VectorDouble vecRhoDt
Time derivative of density.
Definition: Remodeling.hpp:201
MatrixDouble matGradientOfDeformation
Gradient of deformation.
Definition: Remodeling.hpp:198
double getCFromDensity(const double &rho)
Definition: Remodeling.cpp:218
double pSi_ref
reference free energy
Definition: Remodeling.hpp:154
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
Definition: Remodeling.hpp:130
double n
porosity exponent [-]
Definition: Remodeling.hpp:148
int b
b exponent for bell function
Definition: Remodeling.hpp:153
double getCFromDensityDiff(const double &rho)
Definition: Remodeling.cpp:225
double m
algorithmic exponent [-]
Definition: Remodeling.hpp:147
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
Definition: Remodeling.hpp:131
PetscBool is_atom_testing
for atom tests
Definition: Remodeling.hpp:162
double c
density evolution (growth) velocity [d/m^2]
Definition: Remodeling.hpp:146
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
Definition: Remodeling.hpp:136
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
Definition: Remodeling.hpp:138
PetscBool less_post_proc
reduce file size
Definition: Remodeling.hpp:163
DM dm
Discretization manager.
Definition: Remodeling.hpp:127
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
Definition: Remodeling.hpp:141
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
Definition: Remodeling.hpp:137
double R0
mass conduction coefficient
Definition: Remodeling.hpp:155
double cUrrent_psi
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:158
boost::shared_ptr< NonlinearElasticElement > elasticPtr
Definition: Remodeling.hpp:140
double cUrrent_mass
current free energy for evaluating equilibrium state
Definition: Remodeling.hpp:160
MoFEM::Interface & mField
Definition: Remodeling.hpp:206
MoFEMErrorCode addElements()
Set and add finite elements.
MoFEMErrorCode solveDM()
Solve problem set up in DM.
MoFEMErrorCode buildDM()
Set problem and DM.
MoFEMErrorCode getParameters()
Get parameters form line command or config file.
MoFEMErrorCode addFields()
Set and add entities to approximation fields.
MoFEMErrorCode addMomentumFluxes()
Finite elements to calculate tractions.
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, std::string mesh_node_positions="MESH_NODE_POSITIONS")
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
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
Deprecated interface functions.
Basic algebra on fields.
Definition: FieldBlas.hpp:31
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
keeps basic data about problem
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:27
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:166
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
structure grouping operators and data used for calculation of nonlinear elastic element
Force scale operator for reading two columns.