In the following code, we solve the problem of identification of forces which result in given displacements on some surface inside the body.
In the experiment, a body is covered on top by transparent gel layer, on the interface between movements of markers (bids) is observed. Strain in the body is generated by cells living on the top layer. Here we looking for forces generated by those cells.
Following implementing when needed could be easily extended to curved surfaces arbitrary located in the body.
This problem can be mathematically described by minimisation problem with the constraints, as follows
\[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +
\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} \\
\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
n\cdot\sigma(u) = \rho \quad \textrm{on}\,S_\rho \\
\end{array}
\right.
\]
Applying standard procedure, i.e. Lagrange multiplier method, we get
\[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho,\Upsilon) =
\frac{1}{2} (u,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} +
(\Upsilon,\textrm{div}[\sigma])_V \\ n\cdot\sigma = \rho \quad
\textrm{on}\,S_\rho \end{array} \right. \]
\[ J(u,\rho,\Upsilon) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,n\cdot\sigma)_{S_\rho}
\]
and using second constraint, i.e. static constrain
\[
J(u,\rho,\Upsilon) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,\rho)_{S_\rho}.
\]
Note that force is enforced in week sense.
\[
(\textrm{grad}[\Upsilon],\sigma)_V =
(\textrm{grad}[\Upsilon],\mathcal{A} \textrm{grad}[u] )_V =
(\mathcal{A}^{*}\textrm{grad}[\Upsilon],\textrm{grad}[u] )_V =
(\Sigma,\textrm{grad}[u] )_V
\]
\[
J(u,\rho,\Upsilon) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+(\Sigma,\textrm{grad}[u])_V - (\Upsilon,\rho)_{S_\rho}
\]
\[
\left\{
\begin{array}{ll}
(Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V &= 0 \\
(\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} &= 0
\\
(\epsilon_\rho \rho,\delta \rho)_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} &= 0
\end{array}
\right.
\]
and applying finite element discretization, above equations are expressed by linear algebraic equations suitable for computer calculations
\[
\left[
\begin{array}{ccc}
S & K_{u \Upsilon} & 0 \\
K_{\Upsilon u} & 0 & -B^\textrm{T} \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
S u_d \\ 0 \\ 0
\end{array}
\right]
\]
\[
\left[
\begin{array}{ccc}
K & 0 & -B^\textrm{T} \\
S & K & 0 \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
0 \\ S u_d \\ 0
\end{array}
\right]
\]
\[
S=\epsilon_u^{-1} I\\
\epsilon_\rho^\textrm{abs} = \epsilon_\rho^\text{ref}\epsilon_u^{-1}
\]
Displacement field is \(\mathbf{u}\) and \(\Upsilon\), force field on top surface is \(\boldsymbol\rho\).
We not yet explored nature of cell forces. Since the cells are attached to the body surface and anything else, the body as results those forces can not be subjected to rigid body motion.
We assume that forces are generated by 2d objects, as results only tangential forces can be produced by those forces if the surface is planar. For non-planar surfaces addition force normal to the surface is present, although the magnitude of this force could be calculated purely from geometrical considerations, thus additional physical equations are not needed.
Assuming that straight and very short pre-stressed fibres generate cell forces; a force field is curl-free. Utilising that forces can be expressed by potential field as follows
where \(\Phi\) is force potential field.
For this problem \(H^(\Omega)\) function space is required for \(u\), \(\Upsilon\) and \(\rho\). Note that \(\rho\) is given by derivative of potential field \(\Phi\) which derivatives gives potential forces.
In this variant generalisation of the local model to a weakly enforced force curl-free cell force field is considered. This generalisation is a consequence of the observation that cell has some small but finite size. Moreover, a cell has a complex pre-stressed structure which can transfer shear forces.
\[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +
\frac{\epsilon_\rho}{2}(\rho,\rho-l^2\textrm{curl}^s
[\textrm{curl}^s[\rho]])_{S_\rho} \\
\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
n\cdot\sigma = \rho
\end{array}
\right.
\]
\[
\left\{
\begin{array}{l}
\textrm{min}\,J(u,\rho) =
\frac{1}{2} (u-u_d,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
+\frac{\epsilon_l^{-1}}{2}(\textrm{curl}^s[\rho],\textrm{curl}^s[\rho])_{S_\rho}
\\
\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
n\cdot\sigma = \rho
\end{array}
\right.
\]
where parameter \(\epsilon_l\) controls curl-free term, if this parameter approach zero, this formulation converge to local variant presented above.
\[
\left\{
\begin{array}{l}
(Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V= 0 \\
(\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} = 0 \\
(\epsilon_\rho \rho,\delta
\rho)_{S_\rho}+\epsilon_l^{-1}(\textrm{curl}^s[\rho],\textrm{curl}^s[\delta
\rho])_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} = 0 \end{array} \right. \]
\[ \left[ \begin{array}{ccc}
S & K_{u \Upsilon} & 0 \\
K_{\Upsilon u} & 0 & -B^\textrm{T} \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
S u_d \\ 0 \\ 0
\end{array}
\right]
\]
\[
\left[
\begin{array}{ccc}
K & 0 & -B^\textrm{T} \\
S & K & 0 \\
0 & -B & D
\end{array}
\right]
\left[
\begin{array}{c}
u \\ \Upsilon \\ \rho
\end{array}
\right]
=
\left[
\begin{array}{c}
0 \\ S u_d \\ 0
\end{array}
\right]
\]
For this problem \(H^(\Omega)\) function space is required for \(u\), \(\Upsilon\). Note that force field need to be in \(H(\textrm{curl},S_\rho)\) space.
If only one layer is specified, another thin polymer layer could be created automatically using prism elements,
static char help[] =
"-my_block_config set block data\n"
"\n";
};
}
using namespace boost::numeric;
#include <boost/program_options.hpp>
namespace po = boost::program_options;
int main(
int argc,
char *argv[]) {
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
PetscBool flg_block_config, flg_file;
char block_config_file[255];
PetscBool flg_order_force;
PetscInt order_force = 2;
PetscBool flg_eps_u, flg_eps_rho, flg_eps_l;
double eps_u = 1e-6;
double eps_rho = 1e-3;
double eps_l = 0;
PetscBool is_curl = PETSC_TRUE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
"-my_order_force",
"default approximation order for force approximation", "", order_force,
&order_force, &flg_order_force);
CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
"", "block_conf.in", block_config_file, 255,
&flg_block_config);
CHKERR PetscOptionsReal(
"-my_eps_rho",
"foce optimality parameter ",
"",
eps_rho, &eps_rho, &flg_eps_rho);
CHKERR PetscOptionsReal(
"-my_eps_u",
"displacement optimality parameter ",
"", eps_u, &eps_u, &flg_eps_u);
CHKERR PetscOptionsReal(
"-my_eps_l",
"displacement optimality parameter ",
"", eps_l, &eps_l, &flg_eps_l);
CHKERR PetscOptionsBool(
"-my_curl",
"use Hcurl space to approximate forces",
"", is_curl, &is_curl, NULL);
ierr = PetscOptionsEnd();
if (flg_file != PETSC_TRUE) {
"*** ERROR -my_file (MESH FILE NEEDED)");
}
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR mmanager_ptr->printDisplacementSet();
CHKERR mmanager_ptr->printForceSet();
CHKERR mmanager_ptr->printMaterialsSet();
bit_level0.set(0);
{
CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
}
std::vector<Range> setOrderToEnts(10);
std::map<int, BlockOptionData> block_data;
if (flg_block_config) {
double read_eps_u, read_eps_rho, read_eps_l;
try {
ifstream ini_file(block_config_file);
if (!ini_file.is_open()) {
SETERRQ(PETSC_COMM_SELF, 1,
"*** -my_block_config does not exist ***");
}
po::variables_map vm;
po::options_description config_file_options;
config_file_options.add_options()(
"eps_u", po::value<double>(&read_eps_u)->default_value(-1))(
"eps_rho", po::value<double>(&read_eps_rho)->default_value(-1))(
"eps_l", po::value<double>(&read_eps_l)->default_value(-1));
std::ostringstream str_order;
str_order << "block_" << it->getMeshsetId() << ".displacement_order";
config_file_options.add_options()(
str_order.str().c_str(),
po::value<int>(&block_data[it->getMeshsetId()].oRder)
std::ostringstream str_cond;
str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
config_file_options.add_options()(
str_cond.str().c_str(),
po::value<double>(&block_data[it->getMeshsetId()].yOung)
->default_value(-1));
std::ostringstream str_capa;
str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
config_file_options.add_options()(
str_capa.str().c_str(),
po::value<double>(&block_data[it->getMeshsetId()].pOisson)
->default_value(-2));
}
po::parsed_options parsed =
parse_config_file(ini_file, config_file_options, true);
store(parsed, vm);
po::notify(vm);
if (block_data[it->getMeshsetId()].oRder == -1)
continue;
if (block_data[it->getMeshsetId()].oRder ==
order)
continue;
PetscPrintf(PETSC_COMM_WORLD, "Set block %d order to %d\n",
it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
true);
CHKERR moab.get_connectivity(block_ents, nodes,
true);
Range ents_to_set_order, ents3d;
CHKERR moab.get_adjacencies(nodes, 3,
false, ents3d,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(ents3d, 2,
false, ents_to_set_order,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(ents3d, 1,
false, ents_to_set_order,
moab::Interface::UNION);
ents_to_set_order = subtract(
ents_to_set_order, ents_to_set_order.subset_by_type(MBQUAD));
ents_to_set_order = subtract(
ents_to_set_order, ents_to_set_order.subset_by_type(MBPRISM));
set_order_ents.merge(ents3d);
set_order_ents.merge(ents_to_set_order);
setOrderToEnts[block_data[it->getMeshsetId()].oRder].merge(
set_order_ents);
}
std::vector<std::string> additional_parameters;
additional_parameters =
collect_unrecognized(parsed.options, po::include_positional);
for (std::vector<std::string>::iterator vit =
additional_parameters.begin();
vit != additional_parameters.end(); vit++) {
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"** WARNING Unrecognized option %s\n",
vit->c_str());
}
} catch (const std::exception &ex) {
std::ostringstream ss;
ss << ex.what() << std::endl;
}
if (read_eps_u > 0) {
eps_u = read_eps_u;
};
if (read_eps_rho > 0) {
eps_rho = read_eps_rho;
}
if (read_eps_l > 0) {
eps_l = read_eps_l;
}
}
PetscPrintf(PETSC_COMM_WORLD, "epsU = %6.4e epsRho = %6.4e\n", eps_u,
eps_rho);
if (is_curl) {
} else {
}
if (mmanager_ptr->checkMeshset(202,
SIDESET)) {
ents_1st_layer, true);
vertex_to_fix, false);
CHKERR mmanager_ptr->getEntitiesByDimension(202,
SIDESET, 1, edges_to_fix,
false);
if (vertex_to_fix.size() != 1 && !vertex_to_fix.empty()) {
"Should be one vertex only, but is %d", vertex_to_fix.size());
}
}
ents_1st_layer.subset_by_type(MBTRI), MBTRI, "RHO");
if (mmanager_ptr->checkMeshset(101,
SIDESET)) {
ents_2nd_layer, true);
}
for (int oo = 2; oo != setOrderToEnts.size(); oo++) {
if (setOrderToEnts[oo].size() > 0) {
}
}
const int through_thickness_order = 2;
{
CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
CHKERR moab.get_adjacencies(ents3d, 2,
false, ents,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(ents3d, 1,
false, ents,
moab::Interface::UNION);
CHKERR moab.get_entities_by_type(0, MBPRISM, prisms);
{
CHKERR moab.get_adjacencies(prisms, 2,
false, quads,
moab::Interface::UNION);
prism_tris = quads.subset_by_type(MBTRI);
quads = subtract(quads, prism_tris);
CHKERR moab.get_adjacencies(quads, 1,
false, quads_edges,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(prism_tris, 1,
false, prism_tris_edges,
moab::Interface::UNION);
quads_edges = subtract(quads_edges, prism_tris_edges);
prisms.merge(quads);
prisms.merge(quads_edges);
}
ents.merge(ents3d);
ents = subtract(ents, set_order_ents);
ents = subtract(ents, prisms);
through_thickness_order);
}
if (is_curl) {
} else {
}
boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(
new Hooke<adouble>());
boost::shared_ptr<Hooke<double>> hooke_double_ptr(
new Hooke<double>());
CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
CHKERR elastic.addElement(
"ELASTIC",
"U");
CHKERR elastic.setOperators(
"U",
"MESH_NODE_POSITIONS",
false,
true);
{
elastic.setOfBlocks[2].tEts);
elastic.setOfBlocks[2].tEts, MBPRISM, "ELASTIC_PRISM");
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
fat_prism_rhs.getOpPtrVector().push_back(
fat_prism_rhs.getOpPtrVector().push_back(
fat_prism_rhs.getOpPtrVector().push_back(
"U", elastic.commonData));
fat_prism_rhs.getOpPtrVector().push_back(
"U", elastic.setOfBlocks[2], elastic.commonData, 2, false, false,
true));
fat_prism_rhs.getOpPtrVector().push_back(
"U", elastic.setOfBlocks[2], elastic.commonData));
fat_prism_lhs.getOpPtrVector().push_back(
fat_prism_lhs.getOpPtrVector().push_back(
fat_prism_lhs.getOpPtrVector().push_back(
"U", elastic.commonData));
fat_prism_lhs.getOpPtrVector().push_back(
"U", elastic.setOfBlocks[2], elastic.commonData, 2, true, false,
true));
fat_prism_lhs.getOpPtrVector().push_back(
"U", "U", elastic.setOfBlocks[2], elastic.commonData));
}
int id = it->getMeshsetId();
if (block_data[id].yOung > 0) {
elastic.setOfBlocks[id].E = block_data[id].yOung;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Block %d set Young modulus %3.4g\n", id,
elastic.setOfBlocks[id].E);
}
if (block_data[id].pOisson >= -1) {
elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Block %d set Poisson ratio %3.4g\n", id,
elastic.setOfBlocks[id].PoissonRatio);
}
}
"UPSILON");
"U");
"UPSILON");
"U");
"DISP_X");
"DISP_Y");
"DISPLACEMENTS_PENALTY");
"BT");
"B");
"D");
DMType dm_name = "MOFEM";
DM dm_control;
{
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
CHKERR DMSetType(dm_control, dm_name);
bit_level0);
CHKERR DMSetFromOptions(dm_control);
}
ublas::matrix<Mat> nested_matrices(2, 2);
ublas::vector<IS> nested_is_rows(2);
ublas::vector<IS> nested_is_cols(2);
for (
int i = 0;
i != 2;
i++) {
nested_is_rows[
i] = PETSC_NULL;
nested_is_cols[
i] = PETSC_NULL;
for (
int j = 0;
j != 2;
j++) {
nested_matrices(
i,
j) = PETSC_NULL;
}
}
ublas::matrix<Mat> sub_nested_matrices(2, 2);
ublas::vector<IS> sub_nested_is_rows(2);
ublas::vector<IS> sub_nested_is_cols(2);
for (
int i = 0;
i != 2;
i++) {
sub_nested_is_rows[
i] = PETSC_NULL;
sub_nested_is_cols[
i] = PETSC_NULL;
for (
int j = 0;
j != 2;
j++) {
sub_nested_matrices(
i,
j) = PETSC_NULL;
}
}
DM dm_sub_volume_control;
{
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_volume_control);
CHKERR DMSetType(dm_sub_volume_control, dm_name);
"SUB_CONTROL_PROB");
CHKERR DMSetUp(dm_sub_volume_control);
boost::shared_ptr<Problem::SubProblemData> sub_data =
CHKERR sub_data->getRowIs(&nested_is_rows[0]);
CHKERR sub_data->getColIs(&nested_is_cols[0]);
nested_matrices(0, 0) = PETSC_NULL;
}
{
DM dm_sub_sub_elastic;
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_sub_elastic);
CHKERR DMSetType(dm_sub_sub_elastic, dm_name);
"ELASTIC_PROB");
CHKERR DMSetUp(dm_sub_sub_elastic);
Mat Kuu;
CHKERR DMCreateMatrix(dm_sub_sub_elastic, &Kuu);
CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Du);
CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Fu);
SCATTER_REVERSE);
boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
dirihlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
dirihlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
dirihlet_bc_ptr.get());
elastic.getLoopFeRhs().snes_f = Fu;
fat_prism_rhs.snes_f = Fu;
&elastic.getLoopFeRhs());
&fat_prism_rhs);
elastic.getLoopFeLhs().snes_B = Kuu;
fat_prism_lhs.snes_B = Kuu;
&elastic.getLoopFeLhs());
&fat_prism_lhs);
dirihlet_bc_ptr.get());
CHKERR VecGhostUpdateBegin(Fu, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(Fu, ADD_VALUES, SCATTER_REVERSE);
boost::shared_ptr<Problem::SubProblemData> sub_data =
CHKERR sub_data->getRowIs(&sub_nested_is_rows[0]);
CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
sub_nested_matrices(0, 0) = Kuu;
IS isUpsilon;
->isCreateFromProblemFieldToOtherProblemField(
"ELASTIC_PROB",
"U",
ROW,
"SUB_CONTROL_PROB",
"UPSILON",
ROW,
PETSC_NULL, &isUpsilon);
sub_nested_is_rows[1] = isUpsilon;
sub_nested_is_cols[1] = isUpsilon;
sub_nested_matrices(1, 1) = Kuu;
PetscObjectReference((PetscObject)Kuu);
PetscObjectReference((PetscObject)isUpsilon);
cerr << "Kuu" << endl;
MatView(Kuu, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
CHKERR DMDestroy(&dm_sub_sub_elastic);
}
{
DM dm_sub_disp_penalty;
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_disp_penalty);
CHKERR DMSetType(dm_sub_disp_penalty, dm_name);
"S_PROB");
CHKERR DMSetUp(dm_sub_disp_penalty);
Mat S;
CHKERR DMCreateMatrix(dm_sub_disp_penalty, &S);
face_element.getOpPtrVector().push_back(
new OpCellS(S, eps_u));
"DISPLACEMENTS_PENALTY", &face_element);
CHKERR MatAssemblyBegin(S, MAT_FLUSH_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FLUSH_ASSEMBLY);
cerr << "S" << endl;
MatView(S, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
boost::shared_ptr<Problem::SubProblemData> sub_data =
sub_nested_matrices(1, 0) = S;
CHKERR DMDestroy(&dm_sub_disp_penalty);
}
{
DM dm_sub_force_penalty;
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force_penalty);
CHKERR DMSetType(dm_sub_force_penalty, dm_name);
CHKERR DMSetUp(dm_sub_force_penalty);
CHKERR DMCreateMatrix(dm_sub_force_penalty, &
D);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
{
face_d_matrix.getOpPtrVector().push_back(
face_d_matrix.getOpPtrVector().push_back(
if (is_curl) {
face_d_matrix.getOpPtrVector().push_back(
face_d_matrix.getOpPtrVector().push_back(
} else {
face_d_matrix.getOpPtrVector().push_back(
face_d_matrix.getOpPtrVector().push_back(
}
&face_d_matrix);
}
CHKERR MatAssemblyBegin(
D, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
D, MAT_FINAL_ASSEMBLY);
if (is_curl == PETSC_FALSE) {
int nb_dofs_to_fix = 0;
int index_to_fix = 0;
if (!vertex_to_fix.empty()) {
boost::shared_ptr<NumeredDofEntity> dof_ptr;
dof_ptr);
if (dof_ptr) {
nb_dofs_to_fix = 1;
index_to_fix = dof_ptr->getPetscGlobalDofIdx();
cerr << *dof_ptr << endl;
}
}
}
CHKERR MatZeroRowsColumns(
D, nb_dofs_to_fix, &index_to_fix,
eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
} else {
std::vector<int> dofs_to_fix;
for (auto p_eit = edges_to_fix.pair_begin();
p_eit != edges_to_fix.pair_end(); ++p_eit) {
auto lo = row_dofs->lower_bound(
FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
auto hi =
row_dofs->upper_bound(FieldEntity::getHiLocalEntityBitNumber(
bit_number, p_eit->second));
for (; lo != hi; ++lo)
dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
}
CHKERR MatZeroRowsColumns(
D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
}
cerr << "D" << endl;
MatView(
D, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
boost::shared_ptr<Problem::SubProblemData> sub_data =
CHKERR sub_data->getRowIs(&nested_is_rows[1]);
CHKERR sub_data->getColIs(&nested_is_cols[1]);
nested_matrices(1, 1) =
D;
CHKERR DMDestroy(&dm_sub_force_penalty);
}
{
DM dm_sub_force;
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
CHKERR DMSetType(dm_sub_force, dm_name);
Mat UB, UPSILONB;
CHKERR DMCreateMatrix(dm_sub_force, &UB);
CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
CHKERR MatZeroEntries(UPSILONB);
{
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
face_b_matrices.getOpPtrVector().push_back(
face_b_matrices.getOpPtrVector().push_back(
if (is_curl) {
face_b_matrices.getOpPtrVector().push_back(
face_b_matrices.getOpPtrVector().push_back(
face_b_matrices.getOpPtrVector().push_back(
new OpCellCurlB(UB,
"U"));
face_b_matrices.getOpPtrVector().push_back(
} else {
face_b_matrices.getOpPtrVector().push_back(
face_b_matrices.getOpPtrVector().push_back(
face_b_matrices.getOpPtrVector().push_back(
}
}
CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
if (is_curl == PETSC_FALSE) {
int nb_dofs_to_fix = 0;
int index_to_fix = 0;
if (!vertex_to_fix.empty()) {
boost::shared_ptr<NumeredDofEntity> dof_ptr;
dof_ptr);
if (dof_ptr) {
nb_dofs_to_fix = 1;
index_to_fix = dof_ptr->getPetscGlobalDofIdx();
cerr << *dof_ptr << endl;
}
}
}
CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
PETSC_NULL);
CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
PETSC_NULL, PETSC_NULL);
} else {
std::vector<int> dofs_to_fix;
for (auto p_eit = edges_to_fix.pair_begin();
p_eit != edges_to_fix.pair_end(); ++p_eit) {
auto lo = row_dofs->lower_bound(
FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
auto hi =
row_dofs->upper_bound(FieldEntity::getHiLocalEntityBitNumber(
bit_number, p_eit->second));
for (; lo != hi; ++lo)
dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
}
CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
PETSC_NULL, PETSC_NULL);
CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
0, PETSC_NULL, PETSC_NULL);
}
Mat UBT;
CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
cerr << "UBT" << endl;
MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
boost::shared_ptr<Problem::SubProblemData> sub_data =
nested_matrices(0, 1) = UBT;
cerr << "UPSILONB" << endl;
MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
nested_matrices(1, 0) = UPSILONB;
CHKERR DMDestroy(&dm_sub_force);
}
Mat SubA;
CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
&sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
&SubA);
nested_matrices(0, 0) = SubA;
CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
cerr << "Nested SubA" << endl;
MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
}
CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
&nested_is_cols[0], &nested_matrices(0, 0), &
A);
CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
cerr << "Nested A" << endl;
MatView(
A, PETSC_VIEWER_STDOUT_WORLD);
}
CHKERR DMCreateGlobalVector(dm_control, &
D);
CHKERR DMCreateGlobalVector(dm_control, &
F);
{
face_element.getOpPtrVector().push_back(
new OpGetDispX(common_data));
face_element.getOpPtrVector().push_back(
new OpGetDispY(common_data));
face_element.getOpPtrVector().push_back(
&face_element);
CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
}
KSP solver;
{
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetDM(solver, dm_control);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetDMActive(solver, PETSC_FALSE);
CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
PC pc;
CHKERR PCSetType(pc, PCFIELDSPLIT);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
KSP *sub_ksp;
CHKERR PCFieldSplitGetSubKSP(pc, &
n, &sub_ksp);
{
PC sub_pc_0;
CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
}
} else {
"Only works with pre-conditioner PCFIELDSPLIT");
}
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
SCATTER_REVERSE);
CHKERR VecView(
D, PETSC_VIEWER_DRAW_WORLD);
std::string wait;
std::cin >> wait;
}
for (
int i = 0;
i != 2;
i++) {
if (sub_nested_is_rows[
i]) {
CHKERR ISDestroy(&sub_nested_is_rows[
i]);
}
if (sub_nested_is_cols[
i]) {
CHKERR ISDestroy(&sub_nested_is_cols[
i]);
}
for (
int j = 0;
j != 2;
j++) {
if (sub_nested_matrices(
i,
j)) {
CHKERR MatDestroy(&sub_nested_matrices(
i,
j));
}
}
}
for (
int i = 0;
i != 2;
i++) {
CHKERR ISDestroy(&nested_is_rows[
i]);
}
CHKERR ISDestroy(&nested_is_cols[
i]);
}
for (
int j = 0;
j != 2;
j++) {
if (nested_matrices(
i,
j)) {
CHKERR MatDestroy(&nested_matrices(
i,
j));
}
}
}
CHKERR DMDestroy(&dm_sub_volume_control);
{
CHKERR post_proc.generateReferenceElementMesh();
CHKERR post_proc.addFieldValuesPostProc(
"U");
CHKERR post_proc.addFieldValuesGradientPostProc(
"U");
m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "U",
post_proc.commonData, &elastic.setOfBlocks));
CHKERR post_proc.writeFile(
"out.h5m");
elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
elastic.getLoopFeEnergy().eNergy = 0;
&elastic.getLoopFeEnergy());
PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
elastic.getLoopFeEnergy().eNergy);
}
{
CHKERR post_proc_face.generateReferenceElementMesh();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_face.getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
if (is_curl) {
post_proc_face.getOpPtrVector().push_back(
post_proc_face.getOpPtrVector().push_back(
} else {
post_proc_face.getOpPtrVector().push_back(
post_proc_face.getOpPtrVector().push_back(
}
post_proc_face.getOpPtrVector().push_back(
post_proc_face.mapGaussPts, common_data));
CHKERR post_proc_face.writeFile(
"out_tractions.h5m");
}
CHKERR DMDestroy(&dm_control);
}
return 0;
}
Implementation of linear elastic material.
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ HCURL
field with continuous tangents
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_STD_EXCEPTION_THROW
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
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 DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
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.
#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
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
Calculate and assemble g vector.
Calculate and assemble Z matrix.
Calculate and assemble Z matrix.
Calculate and assemble B matrix.
Calculate and assemble D matrix.
Calculate and assemble S matrix.
Shave results on mesh tags for post-processing.
Set Dirichlet boundary conditions on displacements.
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
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
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
Calculate inverse of jacobian for face element.
Transform local reference derivatives of shape functions to global derivatives.
keeps basic data about problem
MoFEMErrorCode getDofByNameEntAndEntDofIdx(const int field_bit_number, const EntityHandle ent, const int ent_dof_idx, const RowColData row_or_col, boost::shared_ptr< NumeredDofEntity > &dof_ptr) const
get DOFs from problem
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator performs automatic differentiation.
structure grouping operators and data used for calculation of nonlinear elastic element
Operator post-procesing stresses for Hook isotropic material.