#include <boost/program_options.hpp>
using namespace std;
namespace po = boost::program_options;
#include <BasicFiniteElements.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]
PetscBool flg_config = PETSC_FALSE;
is_atom_testing = PETSC_FALSE;
with_adol_c = PETSC_FALSE;
char my_config_file_name[255];
rHo_ref = 1000;
pSi_ref = 2.75e4;
rHo_max = rHo_ref + 0.5 * rHo_ref;
rHo_min = rHo_ref - 0.5 * rHo_ref;
R0 = 0.0;
cUrrent_psi = 0.0;
cUrrent_mass = 0.0;
nOremodellingBlock = false;
less_post_proc = PETSC_FALSE;
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,
&rHo_ref, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-rho_max",
"reference bone density",
"", rHo_max,
&rHo_max, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-rho_min",
"reference bone density",
"", rHo_min,
&rHo_min, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-psi_ref",
"reference energy density",
"", pSi_ref,
&pSi_ref, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-r0",
"mass source",
"", R0, &R0, PETSC_NULL);
CHKERR PetscOptionsBool(
"-my_is_atom_test",
"is used with testing, exit with error when diverged",
"", is_atom_testing, &is_atom_testing, PETSC_NULL);
CHKERR PetscOptionsBool(
"-less_post_proc",
"is used to reduce output file size", "",
less_post_proc, &less_post_proc, PETSC_NULL);
"-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);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Max density rHo_max[kg/m3]: %4.3g\n", rHo_max);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Min density rHo_min[kg/m3]: %4.3g\n", rHo_min);
}
ierr = PetscOptionsEnd();
}
inline double Remodeling::CommonData::getCFromDensity(
const double &
rho) {
double mid_rho = 0.5 * (rHo_max + rHo_min);
double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
return 1 / (1 + (
b != 0) * pow(var_h, 2 *
b));
}
inline double Remodeling::CommonData::getCFromDensityDiff(
const double &
rho) {
double mid_rho = 0.5 * (rHo_max + rHo_min);
double var_h = (
rho - mid_rho) / (rHo_max - mid_rho);
return (
b != 0) * (-2) *
b * pow(var_h, 2 *
b - 1) /
((rHo_max - mid_rho) * pow(pow(var_h, 2 *
b) + 1, 2));
}
OpGetRhoTimeDirevative::OpGetRhoTimeDirevative(
Remodeling::CommonData &common_data)
"RHO", UserDataOperator::OPROW),
commonData(common_data) {}
OpGetRhoTimeDirevative::doWork(
int side, EntityType
type,
const int nb_gauss_pts = data.getN().size1();
if (!nb_gauss_pts)
const int nb_base_functions = data.getN().size2();
auto base_function = data.getFTensor0N();
commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
commonData.data.vecRhoDt.clear();
}
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),
commonData(common_data) {
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
if (!commonData.data.gradDispPtr) {
"Gradient at integration points of displacement is not calculated");
}
auto grad = getFTensor2FromMat<3, 3>(*commonData.data.gradDispPtr);
const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
if (!commonData.data.rhoPtr) {
"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);
MatrixDouble &matF = commonData.data.matGradientOfDeformation;
matF.resize(9, nb_gauss_pts, false);
const double c = commonData.c;
const double n = commonData.n;
const double m = commonData.m;
const double rho_ref = commonData.rHo_ref;
const double psi_ref = commonData.pSi_ref;
const double lambda = commonData.lambda;
const double mu = commonData.mu;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
F(0, 0) += 1;
F(1, 1) += 1;
F(2, 2) += 1;
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;
++psi;
++S;
++P;
++R;
++F;
}
}
OpCalculateStressTangentWithAdolc::OpCalculateStressTangentWithAdolc(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
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);
MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
dP_dF.resize(81, nb_gauss_pts, false);
MatrixDouble &matF = commonData.data.matGradientOfDeformation;
matF.resize(9, nb_gauss_pts, false);
commonData.data.matP.resize(9, nb_gauss_pts, false);
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
commonData.data.vecPsi.resize(nb_gauss_pts, false);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
F(0, 0) += 1;
F(1, 1) += 1;
F(2, 2) += 1;
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) {
"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)};
int r_h = ::hessian(commonData.tAg, 6, &vecC[0], hessian_ptr);
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;
++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) {
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
double def_VAL[9];
bzero(def_VAL, 9 * sizeof(double));
Tag th_piola2;
CHKERR postProcMesh.tag_get_handle(
"Piola_Stress2", 9, MB_TYPE_DOUBLE,
th_piola2, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
Tag th_psi;
CHKERR postProcMesh.tag_get_handle(
"psi", 1, MB_TYPE_DOUBLE, th_psi,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
Tag th_rho;
CHKERR postProcMesh.tag_get_handle(
"RHO", 1, MB_TYPE_DOUBLE, th_rho,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
Tag th_piola1;
Tag principal;
Tag th_prin_stress_vect1, th_prin_stress_vect2, th_prin_stress_vect3;
Tag th_R;
Tag th_grad_rho;
if (!commonData.less_post_proc) {
CHKERR postProcMesh.tag_get_handle(
"Piola_Stress1", 9, MB_TYPE_DOUBLE,
th_piola1, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
CHKERR postProcMesh.tag_get_handle(
"Principal_stress", 3, MB_TYPE_DOUBLE,
principal, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
CHKERR postProcMesh.tag_get_handle(
"Principal_stress_S1", 3, MB_TYPE_DOUBLE,
th_prin_stress_vect1,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle(
"Principal_stress_S2", 3, MB_TYPE_DOUBLE,
th_prin_stress_vect2,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle(
"Principal_stress_S3", 3, MB_TYPE_DOUBLE,
th_prin_stress_vect3,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle(
"R", 1, MB_TYPE_DOUBLE, th_R,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle(
"grad_rho", 3, MB_TYPE_DOUBLE,
th_grad_rho,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
}
const int nb_gauss_pts = commonData.data.gradDispPtr->size2();
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
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));
if (!commonData.less_post_proc) {
for (int ii = 0; ii != 3; ii++) {
for (int jj = 0; jj != 3; jj++) {
stress(ii, jj) = P(ii, jj);
}
}
CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
&stress(0, 0));
val = R;
CHKERR postProcMesh.tag_set_data(th_R, &mapGaussPts[gg], 1, &val);
const double t2[] = {grad_rho(0), grad_rho(1), grad_rho(2)};
CHKERR postProcMesh.tag_set_data(th_grad_rho, &mapGaussPts[gg], 1, t2);
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);
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;
++grad_rho;
}
}
OpCalculateStressTangent::OpCalculateStressTangent(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {
doVertices = true;
doEdges = false;
doQuads = false;
doTris = false;
doTets = false;
doPrisms = false;
}
OpCalculateStressTangent::doWork(
int side, EntityType
type,
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);
MatrixDouble &dP_dF = commonData.data.matPushedMaterialTangent;
dP_dF.resize(81, nb_gauss_pts, false);
MatrixDouble &matF = commonData.data.matGradientOfDeformation;
matF.resize(9, nb_gauss_pts, false);
commonData.data.matP.resize(9, nb_gauss_pts, false);
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
commonData.data.vecPsi.resize(nb_gauss_pts, false);
const double lambda = commonData.lambda;
const double mu = commonData.mu;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
F(0, 0) += 1;
F(1, 1) += 1;
F(2, 2) += 1;
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;
invF(
J,
i) * (invF(
L,
k) * lambda) + invF(
L,
i) * (invF(
J,
k) * coef);
D2(0, 0, 0, 0) += mu;
D2(0, 1, 0, 1) += mu;
D2(0, 2, 0, 2) += mu;
D2(1, 0, 1, 0) += mu;
D2(1, 1, 1, 1) += mu;
D2(1, 2, 1, 2) += mu;
D2(2, 0, 2, 0) += mu;
D2(2, 1, 2, 1) += mu;
D2(2, 2, 2, 2) += mu;
++det_f;
++invF;
++grad;
++F;
++D2;
++P;
++psi;
}
}
OpAssmbleStressRhs::OpAssmbleStressRhs(Remodeling::CommonData &common_data)
"DISPLACEMENTS", UserDataOperator::OPROW),
commonData(common_data) {}
OpAssmbleStressRhs::doWork(
int side, EntityType
type,
const int nb_dofs = data.getIndices().size();
if (!nb_dofs)
nF.resize(nb_dofs, false);
nF.clear();
const int nb_gauss_pts = data.getN().size1();
const int nb_base_functions = data.getN().size2();
auto diff_base_functions = data.getFTensor1DiffN<3>();
auto P = getFTensor2FromMat<3, 3>(commonData.data.matP);
const double n = commonData.n;
const double rho_ref = commonData.rHo_ref;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[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);
}
OpAssmbleRhoRhs::OpAssmbleRhoRhs(Remodeling::CommonData &common_data)
"RHO", UserDataOperator::OPROW),
commonData(common_data) {}
OpAssmbleRhoRhs::doWork(
int side, EntityType
type,
const int nb_dofs = data.getIndices().size();
if (!nb_dofs)
nF.resize(nb_dofs, false);
nF.clear();
const int nb_gauss_pts = data.getN().size1();
const int nb_base_functions = data.getN().size2();
auto base_functions = data.getFTensor0N();
auto diff_base_functions = data.getFTensor1DiffN<3>();
if (!commonData.data.rhoPtr) {
}
if (!commonData.data.gradRhoPtr) {
}
auto grad_rho = getFTensor1FromMat<3>(*commonData.data.gradRhoPtr);
if (commonData.data.vecRhoDt.size() != nb_gauss_pts) {
commonData.data.vecRhoDt.resize(nb_gauss_pts, false);
commonData.data.vecRhoDt.clear();
}
const double R0 = commonData.R0;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[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;
++R;
}
getFEMethod()->ts_F,
nb_dofs, &*data.getIndices().begin(), &*nF.begin(), ADD_VALUES);
}
OpAssmbleRhoLhs_dRho::OpAssmbleRhoLhs_dRho(Remodeling::CommonData &common_data)
"RHO", "RHO", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = true;
}
OpAssmbleRhoLhs_dRho::doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
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++) {
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[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;
}
&*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);
&*col_data.getIndices().begin(), row_nb_dofs,
&*row_data.getIndices().begin(),
&transLocK_rho_rho(0, 0), ADD_VALUES);
}
}
OpAssmbleRhoLhs_dF::OpAssmbleRhoLhs_dF(Remodeling::CommonData &common_data)
"RHO", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = false;
}
OpAssmbleRhoLhs_dF::doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
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++) {
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[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++) {
&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;
}
&*row_data.getIndices().begin(), col_nb_dofs,
&*col_data.getIndices().begin(), &locK_rho_F(0, 0),
ADD_VALUES);
}
template <bool ADOLC>
OpAssmbleStressLhs_dF<ADOLC>::OpAssmbleStressLhs_dF(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", "DISPLACEMENTS", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = true;
}
template <bool ADOLC>
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);
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 &dP_dF = commonData.data.matPushedMaterialTangent;
double t0;
const double n = commonData.n;
const double rho_ref = commonData.rHo_ref;
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);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[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();
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);
row_ind, &transLocK_P_F(0, 0), ADD_VALUES);
}
}
OpAssmbleStressLhs_dRho::OpAssmbleStressLhs_dRho(
Remodeling::CommonData &common_data)
"DISPLACEMENTS", "RHO", UserDataOperator::OPROWCOL),
commonData(common_data) {
sYmm = false;
}
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)
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++) {
double w = getVolume() * getGaussPts()(3, gg);
if (getHoGaussPtsDetJac().size() > 0) {
w *= getHoGaussPtsDetJac()[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;
}
&*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),
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");
if (!commonData.less_post_proc) {
CHKERR postProc.addFieldValuesGradientPostProc(
"DISPLACEMENTS");
}
{
boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
&list_of_operators = postProc.getOpPtrVector();
list_of_operators.push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
list_of_operators.push_back(new OpCalculateScalarFieldGradient<3>(
"RHO", commonData.data.gradRhoPtr));
"DISPLACEMENTS", commonData.data.gradDispPtr));
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.postProcMesh, postProcElastic.mapGaussPts,
"DISPLACEMENTS", sit->second, postProcElastic.commonData));
}
iNit = true;
}
if (mass_postproc) {
int mass_vec_ghost[] = {0};
CHKERR VecCreateGhost(PETSC_COMM_WORLD, (!mField.get_comm_rank()) ? 1 : 0,
1, 1, mass_vec_ghost, &mass_vec);
CHKERR VecZeroEntries(mass_vec);
CHKERR VecDuplicate(mass_vec, &energ_vec);
Remodeling::Fe energyCalc(mField);
energyCalc.getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
energyCalc.getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>(
"RHO", commonData.data.gradRhoPtr));
energyCalc.getOpPtrVector().push_back(
commonData.data.gradDispPtr));
energyCalc.getOpPtrVector().push_back(
new 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 =
fabs(energ_sum - commonData.cUrrent_psi) / energ_sum;
double equilibrium_mass_rate =
fabs(mass_sum - commonData.cUrrent_mass) / mass_sum;
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);
}
commonData.cUrrent_psi = energ_sum;
commonData.cUrrent_mass = mass_sum;
}
}
int step;
#if PETSC_VERSION_LE(3, 8, 0)
CHKERR TSGetTimeStepNumber(ts, &step);
#else
CHKERR TSGetStepNumber(ts, &step);
#endif
if ((step) % pRT == 0) {
ostringstream sss;
sss << "out_" << step << ".h5m";
CHKERR postProc.postProcMesh.write_file(sss.str().c_str(),
"MOAB",
"PARALLEL=WRITE_PART");
if (analysis_mesh_flg) {
ostringstream ttt;
ttt << "analysis_mesh_" << step << ".h5m";
CHKERR mField.get_moab().write_file(ttt.str().c_str(),
"MOAB",
"PARALLEL=WRITE_PART");
}
if (commonData.nOremodellingBlock) {
&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");
commonData.oRder = 2;
CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"", 2,
&commonData.oRder, PETSC_NULL);
ierr = PetscOptionsEnd();
CHKERR mField.get_moab().get_entities_by_type(0, MBTET, commonData.tEts_all);
string name = it->getName();
if (name.compare(0, 14, "NO_REMODELLING") == 0) {
commonData.nOremodellingBlock = true;
EntityHandle meshset = it->getMeshset();
CHKERR this->mField.get_moab().get_entities_by_type(
meshset, MBTET, commonData.tEts_block, true);
commonData.tEts_all =
subtract(commonData.tEts_all, commonData.tEts_block);
}
}
}
commonData.bitLevel.set(0);
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 3, commonData.bitLevel);
int order = commonData.oRder;
CHKERR mField.add_field(
"DISPLACEMENTS",
H1,
bool add_rho_field = false;
if (!mField.check_field("RHO")) {
CHKERR mField.add_ents_to_field_by_type(commonData.tEts_all, MBTET,
"RHO");
CHKERR mField.synchronise_field_entities(
"RHO");
add_rho_field = true;
CHKERR mField.set_field_order(0, MBVERTEX,
"RHO", 1);
CHKERR mField.set_field_order(0, MBEDGE,
"RHO",
order - 1);
CHKERR mField.set_field_order(0, MBTRI,
"RHO",
order - 1);
CHKERR mField.set_field_order(0, MBTET,
"RHO",
order - 1);
}
CHKERR mField.add_ents_to_field_by_type(0, MBTET,
"DISPLACEMENTS");
CHKERR mField.add_ents_to_field_by_type(0, MBTET,
"MESH_NODE_POSITIONS");
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);
{
CHKERR mField.set_field_order(0, MBEDGE,
"MESH_NODE_POSITIONS", 2);
}
CHKERR mField.set_field_order(0, MBVERTEX,
"MESH_NODE_POSITIONS", 1);
if (add_rho_field) {
CHKERR mField.getInterface<FieldBlas>()->setField(commonData.rHo_ref,
MBVERTEX, "RHO");
}
Projection10NodeCoordsOnField ent_method_material(mField,
"MESH_NODE_POSITIONS");
CHKERR mField.loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material);
}
CHKERR mField.modify_finite_element_add_field_row(
"FE_REMODELLING",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_col(
"FE_REMODELLING",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data(
"FE_REMODELLING",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data(
"FE_REMODELLING",
"RHO");
CHKERR mField.modify_finite_element_add_field_data(
"FE_REMODELLING",
"MESH_NODE_POSITIONS");
CHKERR mField.modify_finite_element_add_field_row(
"FE_REMODELLING",
"RHO");
CHKERR mField.modify_finite_element_add_field_col(
"FE_REMODELLING",
"RHO");
CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_all, MBTET,
"FE_REMODELLING");
CHKERR mField.add_finite_element(
"ELASTIC");
CHKERR mField.modify_finite_element_add_field_row(
"ELASTIC",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_col(
"ELASTIC",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data(
"ELASTIC",
"DISPLACEMENTS");
CHKERR mField.modify_finite_element_add_field_data(
"ELASTIC",
"MESH_NODE_POSITIONS");
if (commonData.nOremodellingBlock) {
CHKERR mField.add_ents_to_finite_element_by_type(commonData.tEts_block,
MBTET, "ELASTIC");
}
commonData.elasticMaterialsPtr =
commonData.elasticPtr = boost::shared_ptr<NonlinearElasticElement>(
CHKERR commonData.elasticMaterialsPtr->setBlocks(
commonData.elasticPtr->setOfBlocks);
CHKERR commonData.elasticPtr->addElement(
"ELASTIC",
"DISPLACEMENTS");
CHKERR commonData.elasticPtr->setOperators(
"DISPLACEMENTS", "MESH_NODE_POSITIONS", false, true);
CHKERR mField.build_finite_elements();
CHKERR mField.build_adjacencies(commonData.bitLevel);
commonData.data.rhoPtr = boost::shared_ptr<VectorDouble>(
new VectorDouble());
commonData.data.gradRhoPtr =
commonData.data.gradDispPtr =
boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
&list_of_operators_rhs = commonData.feRhs->getOpPtrVector();
list_of_operators_rhs.push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
list_of_operators_rhs.push_back(
new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
list_of_operators_rhs.push_back(new OpGetRhoTimeDirevative(commonData));
"DISPLACEMENTS", commonData.data.gradDispPtr));
list_of_operators_rhs.push_back(new OpAssmbleStressRhs(commonData));
list_of_operators_rhs.push_back(new OpAssmbleRhoRhs(commonData));
boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
&list_of_operators_lhs = commonData.feLhs->getOpPtrVector();
list_of_operators_lhs.push_back(
new OpCalculateScalarFieldValues("RHO", commonData.data.rhoPtr));
list_of_operators_rhs.push_back(
new OpCalculateScalarFieldGradient<3>("RHO", commonData.data.gradRhoPtr));
"DISPLACEMENTS", commonData.data.gradDispPtr));
if (commonData.with_adol_c) {
list_of_operators_lhs.push_back(
new OpCalculateStressTangentWithAdolc(commonData));
list_of_operators_lhs.push_back(
new OpAssmbleStressLhs_dF<true>(commonData));
} else {
list_of_operators_lhs.push_back(new OpCalculateStressTangent(commonData));
list_of_operators_lhs.push_back(
new OpAssmbleStressLhs_dF<false>(commonData));
}
list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dRho(commonData));
list_of_operators_lhs.push_back(new OpAssmbleRhoLhs_dF(commonData));
list_of_operators_lhs.push_back(new OpAssmbleStressLhs_dRho(commonData));
}
CHKERR MetaNeummanForces::addNeumannBCElements(mField,
"DISPLACEMENTS");
CHKERR MetaNeummanForces::setMomentumFluxOperators(
mField, commonData.neumannForces, PETSC_NULL, "DISPLACEMENTS");
PETSC_NULL, "DISPLACEMENTS");
"DISPLACEMENTS");
for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
commonData.neumannForces.begin();
mit != commonData.neumannForces.end(); mit++) {
mit->second->methodsOp.push_back(
}
for (boost::ptr_map<string, NodalForce>::iterator mit =
commonData.nodalForces.begin();
mit != commonData.nodalForces.end(); mit++) {
mit->second->methodsOp.push_back(
}
for (boost::ptr_map<string, EdgeForce>::iterator mit =
commonData.edgeForces.begin();
mit != commonData.edgeForces.end(); mit++) {
mit->second->methodsOp.push_back(
}
}
commonData.dm_name = "DM_REMODELING";
CHKERR DMCreate(PETSC_COMM_WORLD, &commonData.dm);
CHKERR DMSetType(commonData.dm, commonData.dm_name);
commonData.bitLevel);
CHKERR DMSetFromOptions(commonData.dm);
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
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 mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
problem_ptr->
getName(),
ROW,
"DISPLACEMENTS", 0, 3, &is_disp);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
problem_ptr->
getName(),
ROW,
"RHO", 0, 1, &is_rho);
CHKERR PCFieldSplitSetIS(pc, NULL, is_disp);
CHKERR PCFieldSplitSetIS(pc, NULL, is_rho);
}
TsCtx *ts_ctx;
{
ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
}
#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);
if (commonData.is_atom_testing) {
if (commonData.cUrrent_psi < 1.67 || commonData.cUrrent_psi > 1.68)
}
}
int row_side, EntityType row_type,
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);
if (getHoGaussPtsDetJac().size() > 0) {
vol *= getHoGaussPtsDetJac()[gg];
}
energy += vol * psi;
++psi;
}
CHKERR VecSetValue(energVec, 0, energy, ADD_VALUES);
CHKERR VecSetValue(massVec, 0, mass, ADD_VALUES);
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define SYMMETRIC_TENSOR2_VEC_PTR(VEC)
#define SYMMETRIC_TENSOR2_MAT_PTR(MAT)
#define TENSOR2_MAT_PTR(MAT)
#define SYMMETRIC_TENSOR4_MAT_PTR(MAT)
#define TENSOR4_MAT_PTR(MAT)
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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 ...
@ 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 ...
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
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.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< double, DoubleAllocator > VectorDouble
MatrixBoundedArray< double, 9 > MatrixDouble3by3
VectorBoundedArray< double, 3 > VectorDouble3
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
implementation of Data Operators for Forces and Sources
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.
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DEPRECATED PetscErrorCode f_TSMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
DeprecatedCoreInterface Interface
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
const double D
diffusivity
DataForcesAndSourcesCore::EntData EntData
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Deprecated interface functions.
keeps basic data about problem
std::string getName() const
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MonitorPostProc(MoFEM::Interface &m_field, std::map< int, NonlinearElasticElement::BlockData > &set_of_blocks, NonlinearElasticElement::MyVolumeFE &fe_elastic_energy, ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, DataAtIntegrationPts &common_data, BlockData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Force scale operator for reading two columns.