#include <boost/program_options.hpp>
namespace po = boost::program_options;
#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]
PetscBool flg_config = PETSC_FALSE;
char my_config_file_name[255];
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling parameters",
"none");
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;
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;
}
}
CHKERR PetscOptionsScalar(
"-young_modulus",
"get young modulus",
"",
CHKERR PetscOptionsScalar(
"-poisson_ratio",
"get poisson_ratio",
"",
CHKERR PetscOptionsScalar(
"-c",
"density evolution (growth) velocity [d/m^2]",
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,
CHKERR PetscOptionsScalar(
"-rho_max",
"reference bone density",
"",
rHo_max,
CHKERR PetscOptionsScalar(
"-rho_min",
"reference bone density",
"",
rHo_min,
CHKERR PetscOptionsScalar(
"-psi_ref",
"reference energy density",
"",
pSi_ref,
CHKERR PetscOptionsScalar(
"-r0",
"mass source",
"",
R0, &
R0, PETSC_NULL);
CHKERR PetscOptionsBool(
"-my_is_atom_test",
"is used with testing, exit with error when diverged",
CHKERR PetscOptionsBool(
"-less_post_proc",
"is used to reduce output file size", "",
"-with_adolc",
"calculate the material tangent with automatic differentiation", "",
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);
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();
}
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));
}
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),
DataForcesAndSourcesCore::EntData &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();
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;
}
for (; bb != nb_base_functions; bb++)
++base_function;
++drho_dt;
}
CHKERR VecRestoreArray(getFEMethod()->ts_u_t, &array);
}
"DISPLACEMENTS", UserDataOperator::OPROW),
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
DataForcesAndSourcesCore::EntData &data) {
if (type != MBVERTEX)
"Gradient at integration points of displacement is not calculated");
}
"Density at integration points is not calculated");
}
matS.resize(nb_gauss_pts, 6, false);
matF.resize(9, nb_gauss_pts, false);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
CHKERR determinantTensor3by3(
F, det_f);
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 *
lambda * log_det * log_det;
const double coef =
lambda * log_det -
mu;
S(
i,
I) = invF(
i,
J) ^ P(
J,
I);
R =
c * kp * (pow(
rho / rho_ref,
n -
m) * psi - psi_ref);
++det_f;
++invF;
++grad;
++psi;
++S;
++P;
}
}
Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
int side,
EntityType type, DataForcesAndSourcesCore::EntData &data) {
if (type != MBVERTEX)
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);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
CHKERR recordFreeEnergy_dC<FTensor::Tensor0<FTensor::PackPtr<double *, 1>>,
if (r_g < 0) {
"ADOL-C function evaluation with error");
}
S(0, 0) *= 2;
S(1, 1) *= 2;
S(2, 2) *= 2;
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)};
if (r_h < 0) {
"ADOL-C function evaluation with error");
}
for (int ii = 0; ii < 6; ii++) {
for (int jj = 0; jj < ii; jj++) {
dS_dC(jj, is + ii) = dS_dC(ii, is + jj);
}
}
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;
++grad;
++S;
++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) {
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
DataForcesAndSourcesCore::EntData &data) {
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);
}
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double val = psi;
CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &val);
CHKERR postProcMesh.tag_set_data(th_rho, &mapGaussPts[gg], 1, &val);
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));
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);
noalias(eigen_vectors) = stress;
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;
++grad_rho;
}
}
Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
DataForcesAndSourcesCore::EntData &data) {
if (type != MBVERTEX)
dP_dF.resize(81, nb_gauss_pts, false);
matF.resize(9, nb_gauss_pts, false);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
CHKERR determinantTensor3by3(
F, det_f);
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 *
lambda * log_det * log_det;
const double var =
lambda * log_det -
mu;
const double coef =
mu -
lambda * log_det;
invF(
J,
i) * (invF(
L,
k) *
lambda) + invF(
L,
i) * (invF(
J,
k) * coef);
++det_f;
++invF;
++grad;
++D2;
++P;
++psi;
}
}
"DISPLACEMENTS", UserDataOperator::OPROW),
DataForcesAndSourcesCore::EntData &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>();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
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];
t1(
i) +=
a * P(
i,
I) * diff_base_functions(
I);
++diff_base_functions;
}
for (; bb != nb_base_functions; bb++) {
++diff_base_functions;
}
++P;
}
getFEMethod()->ts_F,
nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
}
"RHO", UserDataOperator::OPROW),
DataForcesAndSourcesCore::EntData &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>();
}
}
}
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double w = getVolume() * getGaussPts()(3, gg);
int bb = 0;
for (; bb != nb_dofs; bb++) {
nF[bb] += w * base_functions * drho_dt;
nF[bb] -= w * base_functions *
R;
nF[bb] -= w * R0 * diff_base_functions(
I) * grad_rho(
I);
++base_functions;
++diff_base_functions;
}
for (; bb != nb_base_functions; bb++) {
++base_functions;
++diff_base_functions;
}
++grad_rho;
++drho_dt;
}
getFEMethod()->ts_F,
nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
}
"RHO", "RHO", UserDataOperator::OPROWCOL),
sYmm = true;
}
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data) {
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)
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;
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++) {
double w = getVolume() * getGaussPts()(3, gg);
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;
}
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);
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),
sYmm = false;
}
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data) {
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)
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();
auto col_diff_base_functions = col_data.getFTensor1DiffN<3>();
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
double w = getVolume() * getGaussPts()(3, gg);
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++) {
&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;
}
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),
sYmm = true;
}
template <bool ADOLC>
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data) {
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);
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();
double t0;
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
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));
a * row_diff_base_functions(
J) * col_diff_base_functions(
L);
t_assemble(
i,
k) += diffDiff(
J,
L) * D2(
i,
J,
k,
L);
if (ADOLC) {
t0 = diffDiff(
J,
L) * S(
J,
L);
t_assemble(0, 0) += t0;
t_assemble(1, 1) += t0;
t_assemble(2, 2) += t0;
}
++col_diff_base_functions;
}
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_diff_base_functions;
}
++D2;
++S;
}
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),
sYmm = false;
}
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data) {
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)
locK_P_RHO.resize(row_nb_dofs, col_nb_dofs, false);
locK_P_RHO.clear();
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++) {
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++) {
&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;
}
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)
postProcElastic(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");
"-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);
"-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");
CHKERR postProc.addFieldValuesGradientPostProc(
"DISPLACEMENTS");
}
{
auto &list_of_operators = postProc.getOpPtrVector();
list_of_operators.push_back(
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 =
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};
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(
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(&energ_vec);
if (equilibrium_flg) {
double equilibrium_rate =
double equilibrium_mass_rate =
if (equilibrium_rate < rate) {
PETSC_COMM_WORLD,
"Energy equilibrium state is achieved! Difference = %0.6g %%. \n",
equilibrium_rate * 100);
}
if (equilibrium_mass_rate < rate) {
PETSC_COMM_WORLD,
"Mass equilibrium state is achieved! Difference = %0.6g %%. \n",
equilibrium_mass_rate * 100);
}
}
}
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 postProc.postProcMesh.write_file(sss.str().c_str(),
"MOAB",
"PARALLEL=WRITE_PART");
if (analysis_mesh_flg) {
ostringstream ttt;
ttt << "analysis_mesh_" << step << ".h5m";
"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,
ierr = PetscOptionsEnd();
string name = it->getName();
if (name.compare(0, 14, "NO_REMODELLING") == 0) {
}
}
}
bool add_rho_field = false;
add_rho_field = true;
}
{
}
if (add_rho_field) {
MBVERTEX, "RHO");
}
"MESH_NODE_POSITIONS");
}
"DISPLACEMENTS");
"DISPLACEMENTS");
"DISPLACEMENTS");
"MESH_NODE_POSITIONS");
"FE_REMODELLING");
"DISPLACEMENTS");
"MESH_NODE_POSITIONS");
MBTET, "ELASTIC");
}
"DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
false, false);
list_of_operators_rhs.push_back(
list_of_operators_rhs.push_back(
list_of_operators_rhs.push_back(
new OpGetRhoTimeDirevative(
commonData));
list_of_operators_rhs.push_back(
new OpAssmbleStressRhs(
commonData));
list_of_operators_rhs.push_back(
new OpAssmbleRhoRhs(
commonData));
false, false);
list_of_operators_lhs.push_back(
list_of_operators_rhs.push_back(
list_of_operators_lhs.push_back(
new OpCalculateStressTangentWithAdolc(
commonData));
list_of_operators_lhs.push_back(
} else {
list_of_operators_lhs.push_back(
new OpCalculateStressTangent(
commonData));
list_of_operators_lhs.push_back(
}
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));
}
CHKERR MetaNeummanForces::addNeumannBCElements(
mField,
"DISPLACEMENTS");
CHKERR MetaNeummanForces::setMomentumFluxOperators(
PETSC_NULL, "DISPLACEMENTS");
"DISPLACEMENTS");
for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
mit->second->methodsOp.push_back(
}
for (boost::ptr_map<string, NodalForce>::iterator mit =
mit->second->methodsOp.push_back(
}
for (boost::ptr_map<string, EdgeForce>::iterator mit =
mit->second->methodsOp.push_back(
}
}
}
}
}
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
SNES snes;
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
IS is_disp, is_rho;
CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
problem_ptr->
getName(),
ROW,
"DISPLACEMENTS", 0, 3, &is_disp);
problem_ptr->
getName(),
ROW,
"RHO", 0, 1, &is_rho);
CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
}
{
}
#if PETSC_VERSION_GE(3, 7, 0)
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
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);
}
}
DataForcesAndSourcesCore::EntData &row_data) {
if (row_type != MBVERTEX)
const int nb_gauss_pts = row_data.getN().size1();
double energy = 0;
double mass = 0;
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double vol = getVolume() * getGaussPts()(3, gg);
energy += vol * psi;
++psi;
}
CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
}
}
static PetscErrorCode ierr
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
#define TENSOR2_MAT_PTR(MAT)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
@ MOFEM_STD_EXCEPTION_THROW
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define TENSOR4_MAT_PTR(MAT)
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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 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 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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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
const double c
speed of light (cm/ns)
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)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
constexpr IntegrationType I
double young_modulus
Young modulus.
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)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpAssmbleRhoRhs(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_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)
OpAssmbleStressRhs(Remodeling::CommonData &common_data)
OpCalculateStress(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateStressTangent(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateStressTangentWithAdolc(Remodeling::CommonData &common_data)
OpGetRhoTimeDirevative(Remodeling::CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, Remodeling::CommonData &common_data)
VectorDouble vecDetF
determinant of F
VectorDouble vecPsi
Elastic energy density.
VectorDouble vecR
Mass sorce.
boost::shared_ptr< MatrixDouble > gradDispPtr
Ptr to gradient of displacements matrix container.
boost::shared_ptr< VectorDouble > rhoPtr
Ptr to density matrix container.
MatrixDouble matPushedMaterialTangent
MatrixDouble matP
1st Piola stress
MatrixDouble matS
2nd Piola stress
boost::shared_ptr< MatrixDouble > gradRhoPtr
Gradient of density field.
MatrixDouble matMaterialTangent
MatrixDouble matInvF
inverse of deformation gradient
VectorDouble vecRhoDt
Time derivative of density.
MatrixDouble matGradientOfDeformation
Gradient of deformation.
double getCFromDensity(const double &rho)
double pSi_ref
reference free energy
boost::shared_ptr< Fe > feLhs
FE to make left hand side.
double n
porosity exponent [-]
int b
b exponent for bell function
double getCFromDensityDiff(const double &rho)
double m
algorithmic exponent [-]
DMType dm_name
dm (problem) name
boost::shared_ptr< Fe > feRhs
FE to make right hand side.
PetscBool is_atom_testing
for atom tests
double c
density evolution (growth) velocity [d/m^2]
boost::ptr_map< string, NeummanForcesSurface > neumannForces
Forces on surface.
MoFEMErrorCode getParameters()
boost::ptr_map< string, EdgeForce > edgeForces
Forces on edges.
double rHo_max
max density
PetscBool less_post_proc
reduce file size
DM dm
Discretization manager.
boost::shared_ptr< ElasticMaterials > elasticMaterialsPtr
boost::ptr_map< string, NodalForce > nodalForces
Nodal forces.
double rHo_min
min density
double rHo_ref
reference density
double lambda
Lame parameter.
double R0
mass conduction coefficient
double cUrrent_psi
current free energy for evaluating equilibrium state
boost::shared_ptr< NonlinearElasticElement > elasticPtr
double cUrrent_mass
current free energy for evaluating equilibrium state
MoFEM::Interface & mField
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.
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.
Section manager is used to create indexes and sections.
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.
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
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.