v0.14.0
Searching...
No Matches
cell_forces.cpp
Go to the documentation of this file.
1/**
2 * \file cell_forces.cpp
3 * \brief Calculate cell forces
4 * \example cell_forces.cpp
5
6\tableofcontents
7
8\section cell_potential Implementation with forces potential
9
10In the following code, we solve the problem of identification of forces which
11result in given displacements on some surface inside the body.
12
13In the experiment, a body is covered on top by transparent gel layer, on the
14interface between movements of markers (bids) is observed. Strain in the body is
15generated by cells living on the top layer. Here we looking for forces generated
16by those cells.
17
18Following implementing when needed could be easily extended to curved surfaces
19arbitrary located in the body.
20
21\ref mofem_citation
22
23\htmlonly
24<a href="https://doi.org/10.5281/zenodo.439392"><img
26href="https://doi.org/10.5281/zenodo.439395"><img
28\endhtmlonly
29
30\subsection cell_local Local curl-free formulation
31
32This problem can be mathematically described by minimisation problem with the
33constraints, as follows
34\f[
35\left\{
36\begin{array}{l}
37\textrm{min}\,J(u,\rho) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +
38\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} \\
39\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
40n\cdot\sigma(u) = \rho \quad \textrm{on}\,S_\rho \\
41\end{array}
42\right.
43\f]
44Applying standard procedure, i.e. Lagrange multiplier method, we get
45\f[
46\left\{
47\begin{array}{l}
48\textrm{min}\,J(u,\rho,\Upsilon) =
49\frac{1}{2} (u,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho} +
50(\Upsilon,\textrm{div}[\sigma])_V \\ n\cdot\sigma = \rho \quad
51\textrm{on}\,S_\rho \end{array} \right. \f] and applying differentiation by
52parts \f[ J(u,\rho,\Upsilon) = \frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
53+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
55\f]
56and using second constraint, i.e. static constrain
57\f[
58J(u,\rho,\Upsilon) =
59\frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
60+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
62\f]
63Note that force is enforced in week sense.
64
65Notting that
66\f[
67\sigma = \mathcal{A} \left( \textrm{grad}[u] \right)^{s}
68\f]
69and using symmetry of operator \f$\mathcal{A}\f$,
70\f[
75\f]
76we finally get
77\f[
78J(u,\rho,\Upsilon) =
79\frac{1}{2} (u-u_d,S(u-u_d))_{S_u}
80+\frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
82\f]
83
84Calculating derivatives, we can get Euler equations in following form
85\f[
86\left\{
87\begin{array}{ll}
88(Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V &= 0 \\
90\\
91(\epsilon_\rho \rho,\delta \rho)_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} &= 0
92\end{array}
93\right.
94\f]
95and applying finite element discretization, above equations are expressed
96by linear algebraic equations suitable for computer calculations
97\f[
98\left[
99\begin{array}{ccc}
100S & K_{u \Upsilon} & 0 \\
101K_{\Upsilon u} & 0 & -B^\textrm{T} \\
1020 & -B & D
103\end{array}
104\right]
105\left[
106\begin{array}{c}
107u \\ \Upsilon \\ \rho
108\end{array}
109\right]
110=
111\left[
112\begin{array}{c}
113S u_d \\ 0 \\ 0
114\end{array}
115\right]
116\f]
117and finally swapping first two rows, we finally get
118\f[
119\left[
120\begin{array}{ccc}
121K & 0 & -B^\textrm{T} \\
122S & K & 0 \\
1230 & -B & D
124\end{array}
125\right]
126\left[
127\begin{array}{c}
128u \\ \Upsilon \\ \rho
129\end{array}
130\right]
131=
132\left[
133\begin{array}{c}
1340 \\ S u_d \\ 0
135\end{array}
136\right]
137\f]
138where
139\f[
140S=\epsilon_u^{-1} I\\
141\epsilon_\rho^\textrm{abs} = \epsilon_\rho^\text{ref}\epsilon_u^{-1}
142\f]
143
144Displacement field is \f$\mathbf{u}\f$ and \f$\Upsilon\f$, force field on top
145surface is \f$\boldsymbol\rho\f$.
146
147We not yet explored nature of cell forces. Since the cells are attached to the
148body surface and anything else, the body as results those forces can not be
149subjected to rigid body motion.
150
151We assume that forces are generated by 2d objects, as results only tangential
152forces can be produced by those forces if the surface is planar. For non-planar
153surfaces addition force normal to the surface is present, although the magnitude
154of this force could be calculated purely from geometrical considerations, thus
155additional physical equations are not needed.
156
157Assuming that straight and very short pre-stressed fibres generate cell forces;
158a force field is curl-free. Utilising that forces can be expressed by potential
159field as follows
160\f[
161\rho = \frac{\partial \Phi_\rho}{\partial \mathbf{x}}
162\f]
163where \f$\Phi\f$ is force potential field.
164
165\subsubsection cell_local_approx Approximation
166
167For this problem \f$H^(\Omega)\f$ function space is required for \f$u\f$,
168\f$\Upsilon\f$ and \f$\rho\f$. Note that \f$\rho\f$ is given by derivative of
169potential field \f$\Phi\f$ which derivatives gives potential forces.
170
171\subsection cell_nonlocal Nonlocal, i.e. weak curl-free
172
173In this variant generalisation of the local model to a weakly enforced force
174curl-free cell force field is considered. This generalisation is a consequence
175of the observation that cell has some small but finite size. Moreover, a cell
176has a complex pre-stressed structure which can transfer shear forces.
177
178We define model like before, however this time we add additional term,
179\f[
180 \left\{
181\begin{array}{l}
182\textrm{min}\,J(u,\rho) =
183\frac{1}{2} (u-u_d,S(u-u_d))_{S_u} +
184\frac{\epsilon_\rho}{2}(\rho,\rho-l^2\textrm{curl}^s
185[\textrm{curl}^s[\rho]])_{S_\rho} \\
186\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
187n\cdot\sigma = \rho
188\end{array}
189\right.
190\f]
191and reformulating this we have
192\f[
193\left\{
194\begin{array}{l}
195\textrm{min}\,J(u,\rho) =
196\frac{1}{2} (u-u_d,S(u-u_d))_{S_u} + \frac{\epsilon_\rho}{2}(\rho,\rho)_{S_\rho}
197+\frac{\epsilon_l^{-1}}{2}(\textrm{curl}^s[\rho],\textrm{curl}^s[\rho])_{S_\rho}
198\\
199\textrm{div}[\sigma(u)] = 0 \quad \textrm{in}\,V \\
200n\cdot\sigma = \rho
201\end{array}
202\right.
203\f]
204where parameter \f$\epsilon_l\f$ controls curl-free term, if this parameter
205approach zero, this formulation converge to local variant presented above.
206
207Applying the same procedure like before, set of Euler equations is derived
208\f[
209\left\{
210\begin{array}{l}
211(Su,\delta u)_{S_u} + (\Sigma,\textrm{grad}[\delta u])_V= 0 \\
212(\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} = 0 \\
213(\epsilon_\rho \rho,\delta
214\rho)_{S_\rho}+\epsilon_l^{-1}(\textrm{curl}^s[\rho],\textrm{curl}^s[\delta
215\rho])_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} = 0 \end{array} \right. \f] That
216give following system of equations \f[ \left[ \begin{array}{ccc}
217S & K_{u \Upsilon} & 0 \\
218K_{\Upsilon u} & 0 & -B^\textrm{T} \\
2190 & -B & D
220\end{array}
221\right]
222\left[
223\begin{array}{c}
224u \\ \Upsilon \\ \rho
225\end{array}
226\right]
227=
228\left[
229\begin{array}{c}
230S u_d \\ 0 \\ 0
231\end{array}
232\right]
233\f]
234and swapping the first two rows we get
235\f[
236\left[
237\begin{array}{ccc}
238K & 0 & -B^\textrm{T} \\
239S & K & 0 \\
2400 & -B & D
241\end{array}
242\right]
243\left[
244\begin{array}{c}
245u \\ \Upsilon \\ \rho
246\end{array}
247\right]
248=
249\left[
250\begin{array}{c}
2510 \\ S u_d \\ 0
252\end{array}
253\right]
254\f]
255
256\subsubsection cell_nonlocal_approx Approximation
257
258For this problem \f$H^(\Omega)\f$ function space is required for \f$u\f$,
259\f$\Upsilon\f$. Note that force field need to be in
260\f$H(\textrm{curl},S_\rho)\f$ space.
261
262\section cell_solution Field split pre-conditioner
263
264In this implementation, the block structure of the matrix is utilised. The \e
265PCFIELDSPLIT
266(see
267<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html>)
268is utilised. Matrix is stored using nested matrix format \e MATNEST
269(see
270)<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html>).
271
272The sub-netsetd matrix is created containing two upper-left blocks,
273\f[
274X = \left[
275\begin{array}{cc}
276K & 0 \\
277S & K
278\end{array}
279\right]
280\f]
281this netsed matrix is blocked in
282\f[
283A= \left[
284\begin{array}{cc}
285X & V \\
286H & D
287\end{array}
288\right]
289\f]
290where
291\f[
292V= \left[
293\begin{array}{c}
294-B^\textrm{T} \\
2950
296\end{array}
297\right]
298\f]
299and
300\f[
301H= \left[
302\begin{array}{cc}
3030 & -B^\textrm{T}
304\end{array}
305\right]
306\f]
307
308The system associated with X matrix is solved using \e PCFIELDSPLIT and
309multiplicative relaxation scheme. In following implementation sub-matrix, K is
310factorised only once. The system related to matrix \f$A\f$ composed form
311matrices \f$X,D,V,H\f$ is solved using Schur complement. Note that here we
312exploit MoFEM capability to create sub-problems in easy way.
313
314\section cell_running_code Running code
315
316Approximations of displacements
317\code
318./map_disp \
319-my_data_x ./examples/data_x.csv \
320-my_data_y ./examples/data_y.csv \
321-my_file ./examples/mesh_cell_force_map.cub \
322-ksp_type cg -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
323-lambda 0.01 -my_order 3 -scale 0.05 -cube_size 1 -my_nparts 4
324\endcode
325
326If only one layer is specified, another thin polymer layer could be created
327automatically using prism elements,
328\code
329./map_disp_prism \
330-my_thinckness 0.01 \
331-my_data_x ./examples/data_x.csv \
332-my_data_y ./examples/data_y.csv \
333-my_file ./examples/mesh_for_prisms.cub \
334-ksp_type cg -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
335-lambda 0.01 -my_order 3 -scale 0.05 -cube_size 1 -my_nparts 4
336\endcode
337
338Config file
339\include users_modules/cell_engineering/examples/block_config.in
340
341Calculating forces
342\code
343mpirun -np 4 ./cell_forces \
344-my_file analysis_mesh.h5m \
345-my_order 1 -my_order_force 2 -my_max_post_proc_ref_level 0 \
346-my_block_config ./examples/block_config.in \
347-ksp_type fgmres -ksp_monitor \
348-fieldsplit_1_ksp_type fgmres \
349-fieldsplit_1_pc_type lu \
350-fieldsplit_1_pc_factor_mat_solver_package mumps \
351-fieldsplit_1_ksp_max_it 100 \
352-fieldsplit_1_ksp_monitor \
353-fieldsplit_0_ksp_type gmres \
354-fieldsplit_0_ksp_max_it 25 \
355-fieldsplit_0_fieldsplit_0_ksp_type preonly \
356-fieldsplit_0_fieldsplit_0_pc_type lu \
357-fieldsplit_0_fieldsplit_0_pc_factor_mat_solver_package mumps \
358-fieldsplit_0_fieldsplit_1_ksp_type preonly \
359-fieldsplit_0_fieldsplit_1_pc_type lu \
360-fieldsplit_0_fieldsplit_1_pc_factor_mat_solver_package mumps \
361-ksp_atol 1e-6 -ksp_rtol 0 -my_eps_u 1e-4 -my_curl 1
362\endcode
363
364As results of above we get
365\image html cell_engineering_forces_example.gif "Example: results" width=600px
366
367\section cell_install Installation with docker
368
369- First install docker as per the instructions here:
370<https://docs.docker.com/installation/#installation>
371- Install cell user module: \e docker \e pull \e likask/cell_engineering
372
373\todo Improve documentation
374
375\todo Generalization to curved surfaces. That should not be difficult adding
376metric tensors and taking into account curvature.
377
378\todo Instead of elastic material we can use GEL model developed in other
379module. That will allow to take into account drying and other rheological
380effects.
381
382\todo For nonlinear material instead of KSP solver we can add SNES solver
383
384\todo Physical equations for cell mechanics could be added directly to minimized
385functional or by constrains.
386
387*/
388
389/* This file is part of MoFEM.
390
391 * MoFEM is free software: you can redistribute it and/or modify it under
393 * Free Software Foundation, either version 3 of the License, or (at your
394 * option) any later version.
395 *
396 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
397 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
398 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
399 * License for more details.
400 *
401 * You should have received a copy of the GNU Lesser General Public
403
405using namespace MoFEM;
406#include <Hooke.hpp>
407#include <CellForces.hpp>
408
409static char help[] = "-my_block_config set block data\n"
410 "\n";
411
413
415 int oRder;
416 double yOung;
417 double pOisson;
418 BlockOptionData() : oRder(-1), yOung(-1), pOisson(-2) {}
419};
420
421} // namespace CellEngineering
422
423using namespace boost::numeric;
424using namespace CellEngineering;
425
426#include <boost/program_options.hpp>
427using namespace std;
428namespace po = boost::program_options;
429
430static int debug = 0;
431
432int main(int argc, char *argv[]) {
433
434 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
435
436 try {
437
438 moab::Core mb_instance;
439 moab::Interface &moab = mb_instance;
440
441 PetscBool flg_block_config, flg_file;
442 char mesh_file_name[255];
443 char block_config_file[255];
444 PetscBool flg_order_force;
445 PetscInt order = 2;
446 PetscInt order_force = 2;
447 PetscBool flg_eps_u, flg_eps_rho, flg_eps_l;
448 double eps_u = 1e-6;
449 double eps_rho = 1e-3;
450 double eps_l = 0;
451 PetscBool is_curl = PETSC_TRUE;
452 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
453 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
454 mesh_file_name, 255, &flg_file);
455 CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
456 order, &order, PETSC_NULL);
457 CHKERR PetscOptionsInt(
458 "-my_order_force",
459 "default approximation order for force approximation", "", order_force,
460 &order_force, &flg_order_force);
461 CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
462 "", "block_conf.in", block_config_file, 255,
463 &flg_block_config);
464 CHKERR PetscOptionsReal("-my_eps_rho", "foce optimality parameter ", "",
465 eps_rho, &eps_rho, &flg_eps_rho);
466 CHKERR PetscOptionsReal("-my_eps_u", "displacement optimality parameter ",
467 "", eps_u, &eps_u, &flg_eps_u);
468 CHKERR PetscOptionsReal("-my_eps_l", "displacement optimality parameter ",
469 "", eps_l, &eps_l, &flg_eps_l);
470 CHKERR PetscOptionsBool("-my_curl", "use Hcurl space to approximate forces",
471 "", is_curl, &is_curl, NULL);
472 ierr = PetscOptionsEnd();
473 CHKERRQ(ierr);
474
475 // Reade parameters from line command
476 if (flg_file != PETSC_TRUE) {
477 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
478 "*** ERROR -my_file (MESH FILE NEEDED)");
479 }
480
481 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
482 if (pcomm == NULL)
483 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
484
485 // Read mesh to MOAB
486 const char *option;
488 "PARALLEL_RESOLVE_SHARED_ENTS;"
489 "PARTITION=PARALLEL_PARTITION;";
490 // option = "";
492
493 // Create MoFEM (Joseph) database
494 MoFEM::Core core(moab);
495 MoFEM::Interface &m_field = core;
496
497 auto mmanager_ptr = m_field.getInterface<MeshsetsManager>();
498 // print bcs
499 CHKERR mmanager_ptr->printDisplacementSet();
500 CHKERR mmanager_ptr->printForceSet();
501 // print block sets with materials
502 CHKERR mmanager_ptr->printMaterialsSet();
503
504 // stl::bitset see for more details
505 BitRefLevel bit_level0;
506 bit_level0.set(0);
507 {
508 Range ents3d;
509 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
510 CHKERR m_field.seed_ref_level(ents3d, bit_level0, false);
511 }
512
513 // set app. order
514
515 std::vector<Range> setOrderToEnts(10);
516
517 // configure blocks by parsing config file
518 // it allow to set approximation order for each block independently
519 Range set_order_ents;
520 std::map<int, BlockOptionData> block_data;
521 if (flg_block_config) {
523 try {
524 ifstream ini_file(block_config_file);
525 if (!ini_file.is_open()) {
526 SETERRQ(PETSC_COMM_SELF, 1,
527 "*** -my_block_config does not exist ***");
528 }
529 // std::cerr << block_config_file << std::endl;
530 po::variables_map vm;
531 po::options_description config_file_options;
537 std::ostringstream str_order;
538 str_order << "block_" << it->getMeshsetId() << ".displacement_order";
540 str_order.str().c_str(),
541 po::value<int>(&block_data[it->getMeshsetId()].oRder)
542 ->default_value(order));
543 std::ostringstream str_cond;
544 str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
546 str_cond.str().c_str(),
547 po::value<double>(&block_data[it->getMeshsetId()].yOung)
548 ->default_value(-1));
549 std::ostringstream str_capa;
550 str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
552 str_capa.str().c_str(),
553 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
554 ->default_value(-2));
555 }
556 po::parsed_options parsed =
557 parse_config_file(ini_file, config_file_options, true);
558 store(parsed, vm);
559 po::notify(vm);
561 if (block_data[it->getMeshsetId()].oRder == -1)
562 continue;
563 if (block_data[it->getMeshsetId()].oRder == order)
564 continue;
565 PetscPrintf(PETSC_COMM_WORLD, "Set block %d order to %d\n",
566 it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
567 Range block_ents;
568 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
569 true);
570 // block_ents = block_ents.subset_by_type(MBTET);
571 Range nodes;
572 CHKERR moab.get_connectivity(block_ents, nodes, true);
573 Range ents_to_set_order, ents3d;
574 CHKERR moab.get_adjacencies(nodes, 3, false, ents3d,
575 moab::Interface::UNION);
576 CHKERR moab.get_adjacencies(ents3d, 2, false, ents_to_set_order,
577 moab::Interface::UNION);
578 CHKERR moab.get_adjacencies(ents3d, 1, false, ents_to_set_order,
579 moab::Interface::UNION);
580 ents_to_set_order = subtract(
582 ents_to_set_order = subtract(
583 ents_to_set_order, ents_to_set_order.subset_by_type(MBPRISM));
584 set_order_ents.merge(ents3d);
585 set_order_ents.merge(ents_to_set_order);
586 setOrderToEnts[block_data[it->getMeshsetId()].oRder].merge(
587 set_order_ents);
588 }
589 CHKERR m_field.synchronise_entities(set_order_ents, 0);
592 collect_unrecognized(parsed.options, po::include_positional);
593 for (std::vector<std::string>::iterator vit =
595 vit != additional_parameters.end(); vit++) {
596 CHKERR PetscPrintf(PETSC_COMM_WORLD,
597 "** WARNING Unrecognized option %s\n",
598 vit->c_str());
599 }
600 } catch (const std::exception &ex) {
601 std::ostringstream ss;
602 ss << ex.what() << std::endl;
603 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
604 }
605 if (read_eps_u > 0) {
607 };
608 if (read_eps_rho > 0) {
610 }
611 if (read_eps_l > 0) {
613 }
614 }
615
616 PetscPrintf(PETSC_COMM_WORLD, "epsU = %6.4e epsRho = %6.4e\n", eps_u,
617 eps_rho);
618
619 // Fields
620 CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 3, MB_TAG_SPARSE,
621 MF_ZERO);
622 CHKERR m_field.add_field("UPSILON", H1, AINSWORTH_LEGENDRE_BASE, 3,
623 MB_TAG_SPARSE, MF_ZERO);
624 if (is_curl) {
625 CHKERR m_field.add_field("RHO", HCURL, AINSWORTH_LEGENDRE_BASE, 1);
626 } else {
627 CHKERR m_field.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1);
628 }
629
630 // add entitities (by tets) to the field
635
637 CHKERR m_field.synchronise_field_entities("UPSILON");
638 CHKERR m_field.synchronise_field_entities("RHO");
639
640 Range vertex_to_fix;
641 Range edges_to_fix;
642 Range ents_1st_layer;
643 // If problem is partitioned meshset culd not exist on this part
644 if (mmanager_ptr->checkMeshset(202, SIDESET)) {
645 CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 2,
646 ents_1st_layer, true);
647 CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 0,
648 vertex_to_fix, false);
649 CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 1, edges_to_fix,
650 false);
651 if (vertex_to_fix.size() != 1 && !vertex_to_fix.empty()) {
652 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
653 "Should be one vertex only, but is %d", vertex_to_fix.size());
654 }
655 }
656 CHKERR m_field.synchronise_entities(ents_1st_layer, 0);
658 ents_1st_layer.subset_by_type(MBTRI), MBTRI, "RHO");
659 Range ents_2nd_layer;
660 // If problem is partitioned meshset culd not exist on this part
661 if (mmanager_ptr->checkMeshset(101, SIDESET)) {
662 CHKERR mmanager_ptr->getEntitiesByDimension(101, SIDESET, 2,
663 ents_2nd_layer, true);
664 }
665 CHKERR m_field.synchronise_entities(ents_2nd_layer, 0);
666
667 for (int oo = 2; oo != setOrderToEnts.size(); oo++) {
668 if (setOrderToEnts[oo].size() > 0) {
669 CHKERR m_field.synchronise_entities(setOrderToEnts[oo], 0);
670 CHKERR m_field.set_field_order(setOrderToEnts[oo], "U", oo);
671 CHKERR m_field.set_field_order(setOrderToEnts[oo], "UPSILON", oo);
672 }
673 }
674
675 const int through_thickness_order = 2;
676 {
677 Range ents3d;
678 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
679 Range ents;
680 CHKERR moab.get_adjacencies(ents3d, 2, false, ents,
681 moab::Interface::UNION);
682 CHKERR moab.get_adjacencies(ents3d, 1, false, ents,
683 moab::Interface::UNION);
684
685 Range prisms;
686 CHKERR moab.get_entities_by_type(0, MBPRISM, prisms);
687 {
690 moab::Interface::UNION);
691 Range prism_tris;
696 moab::Interface::UNION);
697 Range prism_tris_edges;
698 CHKERR moab.get_adjacencies(prism_tris, 1, false, prism_tris_edges,
699 moab::Interface::UNION);
703 }
704
705 ents.merge(ents3d);
706 ents = subtract(ents, set_order_ents);
707 ents = subtract(ents, prisms);
708
709 CHKERR m_field.synchronise_entities(ents, 0);
710 CHKERR m_field.synchronise_entities(prisms, 0);
711
712 CHKERR m_field.set_field_order(ents, "U", order);
713 CHKERR m_field.set_field_order(ents, "UPSILON", order);
714 // approx. order through thickness to 2
715 CHKERR m_field.set_field_order(prisms, "U", through_thickness_order);
716 CHKERR m_field.set_field_order(prisms, "UPSILON",
717 through_thickness_order);
718 }
719 CHKERR m_field.set_field_order(0, MBVERTEX, "U", 1);
720 CHKERR m_field.set_field_order(0, MBVERTEX, "UPSILON", 1);
721
722 if (is_curl) {
723 CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
724 CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
725 } else {
726 CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
727 CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
728 CHKERR m_field.set_field_order(0, MBVERTEX, "RHO", 1);
729 }
730
733 boost::shared_ptr<Hooke<double>> hooke_double_ptr(new Hooke<double>());
734 NonlinearElasticElement elastic(m_field, 2);
737 CHKERR elastic.setOperators("U", "MESH_NODE_POSITIONS", false, true);
739 CellEngineering::FatPrism fat_prism_rhs(m_field);
740 CellEngineering::FatPrism fat_prism_lhs(m_field);
741 {
746 Range ents_2nd_layer;
747 CHKERR mmanager_ptr->getEntitiesByDimension(2, BLOCKSET, 3,
748 elastic.setOfBlocks[2].tEts);
750 elastic.setOfBlocks[2].tEts, MBPRISM, "ELASTIC_PRISM");
751 // right hand side operators
752 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
753 fat_prism_rhs.getOpPtrVector().push_back(
754 new MoFEM::OpCalculateInvJacForFatPrism(inv_jac_ptr));
755 fat_prism_rhs.getOpPtrVector().push_back(
756 new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac_ptr));
757 fat_prism_rhs.getOpPtrVector().push_back(
759 "U", elastic.commonData));
760 fat_prism_rhs.getOpPtrVector().push_back(
762 "U", elastic.setOfBlocks[2], elastic.commonData, 2, false, false,
763 true));
764 fat_prism_rhs.getOpPtrVector().push_back(
766 "U", elastic.setOfBlocks[2], elastic.commonData));
767 // Left hand side operators
768 fat_prism_lhs.getOpPtrVector().push_back(
769 new MoFEM::OpCalculateInvJacForFatPrism(inv_jac_ptr));
770 fat_prism_lhs.getOpPtrVector().push_back(
771 new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac_ptr));
772 fat_prism_lhs.getOpPtrVector().push_back(
774 "U", elastic.commonData));
775 fat_prism_lhs.getOpPtrVector().push_back(
777 "U", elastic.setOfBlocks[2], elastic.commonData, 2, true, false,
778 true));
779 fat_prism_lhs.getOpPtrVector().push_back(
781 "U", "U", elastic.setOfBlocks[2], elastic.commonData));
782 }
783
784 // Update material parameters
786 int id = it->getMeshsetId();
787 if (block_data[id].yOung > 0) {
788 elastic.setOfBlocks[id].E = block_data[id].yOung;
789 CHKERR PetscPrintf(PETSC_COMM_WORLD,
790 "Block %d set Young modulus %3.4g\n", id,
791 elastic.setOfBlocks[id].E);
792 }
793 if (block_data[id].pOisson >= -1) {
794 elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
795 CHKERR PetscPrintf(PETSC_COMM_WORLD,
796 "Block %d set Poisson ratio %3.4g\n", id,
797 elastic.setOfBlocks[id].PoissonRatio);
798 }
799 }
800
801 // build field
802 CHKERR m_field.build_fields();
803
804 // Control elements
811
814 "UPSILON");
816 "U");
818 "UPSILON");
820 "U");
822 "DISP_X");
824 "DISP_Y");
826 "DISPLACEMENTS_PENALTY");
827
828 // Add element to calculate residual on 1st layer
837 "BT");
838
847 "B");
848
849 // Add element to calculate residual on 1st layer
855 "D");
856
857 // build finite elemnts
861
862 // Register MOFEM DM
863 DMType dm_name = "MOFEM";
864 CHKERR DMRegister_MoFEM(dm_name);
865
866 DM dm_control;
867 {
868 // Craete dm_control instance
869 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
870 CHKERR DMSetType(dm_control, dm_name);
871 // set dm_control data structure which created mofem datastructures
872 CHKERR DMMoFEMCreateMoFEM(dm_control, &m_field, "CONTROL_PROB",
873 bit_level0);
874 CHKERR DMSetFromOptions(dm_control);
875 CHKERR DMMoFEMSetSquareProblem(dm_control, PETSC_TRUE);
876 CHKERR DMMoFEMSetIsPartitioned(dm_control, PETSC_TRUE);
877 // add elements to dm_control
885 CHKERR DMSetUp(dm_control);
886 }
887
888 CellEngineering::CommonData common_data;
889
890 ublas::matrix<Mat> nested_matrices(2, 2);
891 ublas::vector<IS> nested_is_rows(2);
892 ublas::vector<IS> nested_is_cols(2);
893 for (int i = 0; i != 2; i++) {
894 nested_is_rows[i] = PETSC_NULL;
895 nested_is_cols[i] = PETSC_NULL;
896 for (int j = 0; j != 2; j++) {
897 nested_matrices(i, j) = PETSC_NULL;
898 }
899 }
900
901 ublas::matrix<Mat> sub_nested_matrices(2, 2);
902 ublas::vector<IS> sub_nested_is_rows(2);
903 ublas::vector<IS> sub_nested_is_cols(2);
904 for (int i = 0; i != 2; i++) {
905 sub_nested_is_rows[i] = PETSC_NULL;
906 sub_nested_is_cols[i] = PETSC_NULL;
907 for (int j = 0; j != 2; j++) {
908 sub_nested_matrices(i, j) = PETSC_NULL;
909 }
910 }
911
912 DM dm_sub_volume_control;
913 {
914 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_volume_control);
915 CHKERR DMSetType(dm_sub_volume_control, dm_name);
916 // set dm_sub_volume_control data structure which created mofem data
917 // structures
918 CHKERR DMMoFEMCreateSubDM(dm_sub_volume_control, dm_control,
919 "SUB_CONTROL_PROB");
920 CHKERR DMMoFEMSetSquareProblem(dm_sub_volume_control, PETSC_TRUE);
925 // add elements to dm_sub_volume_control
926 CHKERR DMSetUp(dm_sub_volume_control);
927
928 const Problem *prb_ptr;
929 CHKERR m_field.get_problem("SUB_CONTROL_PROB", &prb_ptr);
930 boost::shared_ptr<Problem::SubProblemData> sub_data =
931 prb_ptr->getSubData();
932
933 CHKERR sub_data->getRowIs(&nested_is_rows[0]);
934 CHKERR sub_data->getColIs(&nested_is_cols[0]);
935 // That will be filled at the end
936 nested_matrices(0, 0) = PETSC_NULL;
937 }
938
939 {
940 DM dm_sub_sub_elastic;
941
942 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_sub_elastic);
943 CHKERR DMSetType(dm_sub_sub_elastic, dm_name);
944 // set dm_sub_sub_elastic data structure which created mofem data
945 // structures
946 CHKERR DMMoFEMCreateSubDM(dm_sub_sub_elastic, dm_sub_volume_control,
947 "ELASTIC_PROB");
948 CHKERR DMMoFEMSetSquareProblem(dm_sub_sub_elastic, PETSC_TRUE);
953 // add elements to dm_sub_sub_elastic
954 CHKERR DMSetUp(dm_sub_sub_elastic);
955
956 Mat Kuu;
957 Vec Du, Fu;
958 CHKERR DMCreateMatrix(dm_sub_sub_elastic, &Kuu);
959 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Du);
960 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Fu);
961 CHKERR MatZeroEntries(Kuu);
962 CHKERR VecZeroEntries(Du);
963 CHKERR VecZeroEntries(Fu);
964 CHKERR DMoFEMMeshToLocalVector(dm_sub_sub_elastic, Du, INSERT_VALUES,
965 SCATTER_REVERSE);
966
967 // Manage Dirichlet bC
968 boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
969 dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
970 new DirichletDisplacementBc(m_field, "U", Kuu, Du, Fu));
971 dirihlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
972 dirihlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
973 // preproc
974 CHKERR DMoFEMPreProcessFiniteElements(dm_sub_sub_elastic,
975 dirihlet_bc_ptr.get());
976 // internal force vector (to take into account Dirichlet boundary
977 // conditions)
978 elastic.getLoopFeRhs().snes_f = Fu;
979 fat_prism_rhs.snes_f = Fu;
980 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
981 &elastic.getLoopFeRhs());
982 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
983 &fat_prism_rhs);
984 // elastic element matrix
985 elastic.getLoopFeLhs().snes_B = Kuu;
986 fat_prism_lhs.snes_B = Kuu;
987 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
988 &elastic.getLoopFeLhs());
989 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
990 &fat_prism_lhs);
991 // postproc
992 CHKERR DMoFEMPostProcessFiniteElements(dm_sub_sub_elastic,
993 dirihlet_bc_ptr.get());
996 CHKERR VecAssemblyBegin(Fu);
997 CHKERR VecAssemblyEnd(Fu);
998 CHKERR VecScale(Fu, -1);
999 CHKERR VecDestroy(&Du);
1000 CHKERR VecDestroy(&Fu);
1001
1002 const Problem *prb_ptr;
1003 CHKERR m_field.get_problem("ELASTIC_PROB", &prb_ptr);
1004 boost::shared_ptr<Problem::SubProblemData> sub_data =
1005 prb_ptr->getSubData();
1006
1007 CHKERR sub_data->getRowIs(&sub_nested_is_rows[0]);
1008 CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1009 sub_nested_matrices(0, 0) = Kuu;
1010 IS isUpsilon;
1011 CHKERR m_field.getInterface<ISManager>()
1012 ->isCreateFromProblemFieldToOtherProblemField(
1013 "ELASTIC_PROB", "U", ROW, "SUB_CONTROL_PROB", "UPSILON", ROW,
1014 PETSC_NULL, &isUpsilon);
1015 sub_nested_is_rows[1] = isUpsilon;
1016 sub_nested_is_cols[1] = isUpsilon;
1017 sub_nested_matrices(1, 1) = Kuu;
1018 PetscObjectReference((PetscObject)Kuu);
1019 PetscObjectReference((PetscObject)isUpsilon);
1020
1021 // Matrix View
1022 if (debug) {
1023 cerr << "Kuu" << endl;
1024 MatView(Kuu, PETSC_VIEWER_DRAW_WORLD);
1025 std::string wait;
1026 std::cin >> wait;
1027 }
1028
1029 CHKERR DMDestroy(&dm_sub_sub_elastic);
1030 }
1031
1032 {
1033 DM dm_sub_disp_penalty;
1034
1035 // Craete dm_control instance
1036 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_disp_penalty);
1037 CHKERR DMSetType(dm_sub_disp_penalty, dm_name);
1038 // set dm_sub_disp_penalty data structure which created mofem
1039 // datastructures
1040 CHKERR DMMoFEMCreateSubDM(dm_sub_disp_penalty, dm_sub_volume_control,
1041 "S_PROB");
1042 CHKERR DMMoFEMSetSquareProblem(dm_sub_disp_penalty, PETSC_FALSE);
1045 // add elements to dm_sub_disp_penalty
1047 CHKERR DMSetUp(dm_sub_disp_penalty);
1048
1049 Mat S;
1050 CHKERR DMCreateMatrix(dm_sub_disp_penalty, &S);
1051 CHKERR MatZeroEntries(S);
1052
1053 CellEngineering::FaceElement face_element(m_field);
1054 face_element.getOpPtrVector().push_back(new OpCellS(S, eps_u));
1055 CHKERR DMoFEMLoopFiniteElements(dm_sub_disp_penalty,
1056 "DISPLACEMENTS_PENALTY", &face_element);
1057 CHKERR MatAssemblyBegin(S, MAT_FLUSH_ASSEMBLY);
1058 CHKERR MatAssemblyEnd(S, MAT_FLUSH_ASSEMBLY);
1059
1060 // // Matrix View
1061 if (debug) {
1062 cerr << "S" << endl;
1063 MatView(S, PETSC_VIEWER_DRAW_WORLD);
1064 std::string wait;
1065 std::cin >> wait;
1066 }
1067
1068 const Problem *problem_ptr;
1069 CHKERR m_field.get_problem("S_PROB", &problem_ptr);
1070 boost::shared_ptr<Problem::SubProblemData> sub_data =
1071 problem_ptr->getSubData();
1072 // CHKERR sub_data->getRowIs(&sub_nested_is_rows[1]);
1073 // CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1074 sub_nested_matrices(1, 0) = S;
1075
1076 CHKERR DMDestroy(&dm_sub_disp_penalty);
1077 }
1078
1079 // Calculate penalty matrix
1080 {
1081 DM dm_sub_force_penalty;
1082
1083 // Craete dm_control instance
1084 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force_penalty);
1085 CHKERR DMSetType(dm_sub_force_penalty, dm_name);
1086 // set dm_sub_force_penalty data structure which created mofem
1087 // datastructures
1088 CHKERR DMMoFEMCreateSubDM(dm_sub_force_penalty, dm_control, "D_PROB");
1089 CHKERR DMMoFEMSetSquareProblem(dm_sub_force_penalty, PETSC_TRUE);
1092 // add elements to dm_sub_force_penalty
1094 CHKERR DMSetUp(dm_sub_force_penalty);
1095
1096 Mat D;
1097 CHKERR DMCreateMatrix(dm_sub_force_penalty, &D);
1098 CHKERR MatZeroEntries(D);
1099
1100 auto det_ptr = boost::make_shared<VectorDouble>();
1101 auto jac_ptr = boost::make_shared<MatrixDouble>();
1102 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1103
1104 {
1105 CellEngineering::FaceElement face_d_matrix(m_field);
1106
1107 face_d_matrix.getOpPtrVector().push_back(
1108 new OpCalculateHOJacForFace(jac_ptr));
1109 face_d_matrix.getOpPtrVector().push_back(
1110 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
1111
1112 if (is_curl) {
1113 face_d_matrix.getOpPtrVector().push_back(
1114 new OpSetInvJacHcurlFace(inv_jac_ptr));
1115 face_d_matrix.getOpPtrVector().push_back(
1116 new OpCellCurlD(D, eps_rho / eps_u, eps_l * eps_u));
1117 } else {
1118 face_d_matrix.getOpPtrVector().push_back(
1119 new OpSetInvJacH1ForFace(inv_jac_ptr));
1120 face_d_matrix.getOpPtrVector().push_back(
1121 new OpCellPotentialD(D, eps_rho / eps_u));
1122 }
1123 CHKERR DMoFEMLoopFiniteElements(dm_sub_force_penalty, "D",
1124 &face_d_matrix);
1125 }
1126 CHKERR MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY);
1127 CHKERR MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY);
1128
1129 const Problem *problem_ptr;
1130 CHKERR m_field.get_problem("D_PROB", &problem_ptr);
1131
1132 // Zero rows, force field is given by gradients of potential field, so one
1133 // of the values has to be fixed like for rigid body motion.
1134 if (is_curl == PETSC_FALSE) {
1135 int nb_dofs_to_fix = 0;
1136 int index_to_fix = 0;
1137 if (!vertex_to_fix.empty()) {
1138 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1140 m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
1141 dof_ptr);
1142 if (dof_ptr) {
1143 if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1144 nb_dofs_to_fix = 1;
1145 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1146 cerr << *dof_ptr << endl;
1147 }
1148 }
1149 }
1150 CHKERR MatZeroRowsColumns(D, nb_dofs_to_fix, &index_to_fix,
1151 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1152 } else {
1153 std::vector<int> dofs_to_fix;
1154 for (auto p_eit = edges_to_fix.pair_begin();
1155 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1156 auto bit_number = m_field.get_field_bit_number("RHO");
1157 auto row_dofs = problem_ptr->numeredRowDofsPtr;
1158 auto lo = row_dofs->lower_bound(
1159 FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
1160 auto hi =
1161 row_dofs->upper_bound(FieldEntity::getHiLocalEntityBitNumber(
1162 bit_number, p_eit->second));
1163 for (; lo != hi; ++lo)
1164 if ((*lo)->getPart() == m_field.get_comm_rank())
1165 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1166 }
1167 CHKERR MatZeroRowsColumns(D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1168 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1169 }
1170
1171 // Matrix View
1172 if (debug) {
1173 cerr << "D" << endl;
1174 MatView(D, PETSC_VIEWER_DRAW_WORLD);
1175 std::string wait;
1176 std::cin >> wait;
1177 }
1178
1179 boost::shared_ptr<Problem::SubProblemData> sub_data =
1180 problem_ptr->getSubData();
1181 CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1182 CHKERR sub_data->getColIs(&nested_is_cols[1]);
1183 nested_matrices(1, 1) = D;
1184
1185 CHKERR DMDestroy(&dm_sub_force_penalty);
1186 }
1187
1188 {
1189 DM dm_sub_force;
1190
1191 // Craete dm_control instance
1192 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
1193 CHKERR DMSetType(dm_sub_force, dm_name);
1194 // set dm_sub_force data structure which created mofem data structures
1195 CHKERR DMMoFEMCreateSubDM(dm_sub_force, dm_control, "FORCES_PROB");
1196 CHKERR DMMoFEMSetSquareProblem(dm_sub_force, PETSC_FALSE);
1200 // add elements to dm_sub_force
1202 CHKERR DMSetUp(dm_sub_force);
1203
1204 Mat UB, UPSILONB;
1205 CHKERR DMCreateMatrix(dm_sub_force, &UB);
1206 CHKERR MatZeroEntries(UB);
1207 // note be will be transposed in place latter on
1208 CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
1209 CHKERR MatZeroEntries(UPSILONB);
1210 {
1211 CellEngineering::FaceElement face_b_matrices(m_field);
1212 auto det_ptr = boost::make_shared<VectorDouble>();
1213 auto jac_ptr = boost::make_shared<MatrixDouble>();
1214 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1215 face_b_matrices.getOpPtrVector().push_back(
1216 new OpCalculateHOJacForFace(jac_ptr));
1217 face_b_matrices.getOpPtrVector().push_back(
1218 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
1219 if (is_curl) {
1220 face_b_matrices.getOpPtrVector().push_back(
1221 new OpSetInvJacH1ForFace(inv_jac_ptr));
1222 face_b_matrices.getOpPtrVector().push_back(
1223 new OpSetInvJacHcurlFace(inv_jac_ptr));
1224 face_b_matrices.getOpPtrVector().push_back(new OpCellCurlB(UB, "U"));
1225 face_b_matrices.getOpPtrVector().push_back(
1226 new OpCellCurlB(UPSILONB, "UPSILON"));
1227 } else {
1228 face_b_matrices.getOpPtrVector().push_back(
1229 new OpSetInvJacH1ForFace(inv_jac_ptr));
1230 face_b_matrices.getOpPtrVector().push_back(
1231 new OpCellPotentialB(UB, "U"));
1232 face_b_matrices.getOpPtrVector().push_back(
1233 new OpCellPotentialB(UPSILONB, "UPSILON"));
1234 }
1235 CHKERR DMoFEMLoopFiniteElements(dm_sub_force, "B", &face_b_matrices);
1236 }
1237 CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
1238 CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
1239 CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
1240 CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
1241
1242 const Problem *problem_ptr;
1243 CHKERR m_field.get_problem("FORCES_PROB", &problem_ptr);
1244
1245 // Zero rows, force field is given by gradients of potential field, so one
1246 // of the values has to be fixed like for rigid body motion.
1247 if (is_curl == PETSC_FALSE) {
1248 int nb_dofs_to_fix = 0;
1249 int index_to_fix = 0;
1250 if (!vertex_to_fix.empty()) {
1251 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1253 m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
1254 dof_ptr);
1255 if (dof_ptr) {
1256 if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1257 nb_dofs_to_fix = 1;
1258 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1259 cerr << *dof_ptr << endl;
1260 }
1261 }
1262 }
1263 CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
1264 PETSC_NULL);
1265 CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
1266 PETSC_NULL, PETSC_NULL);
1267 } else {
1268 std::vector<int> dofs_to_fix;
1269 for (auto p_eit = edges_to_fix.pair_begin();
1270 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1271 auto bit_number = m_field.get_field_bit_number("RHO");
1272 auto row_dofs = problem_ptr->numeredRowDofsPtr;
1273 auto lo = row_dofs->lower_bound(
1274 FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
1275 auto hi =
1276 row_dofs->upper_bound(FieldEntity::getHiLocalEntityBitNumber(
1277 bit_number, p_eit->second));
1278 for (; lo != hi; ++lo)
1279 if ((*lo)->getPart() == m_field.get_comm_rank())
1280 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1281 }
1282 CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
1283 PETSC_NULL, PETSC_NULL);
1284 CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1285 0, PETSC_NULL, PETSC_NULL);
1286 }
1287
1288 Mat UBT;
1289 CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
1290 CHKERR MatDestroy(&UB);
1291
1292 // Matrix View
1293 if (debug) {
1294 cerr << "UBT" << endl;
1295 MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
1296 std::string wait;
1297 std::cin >> wait;
1298 }
1299
1300 boost::shared_ptr<Problem::SubProblemData> sub_data =
1301 problem_ptr->getSubData();
1302 // CHKERR sub_data->getColIs(&nested_is_rows[0]);
1303 // CHKERR sub_data->getRowIs(&nested_is_cols[1]);
1304 nested_matrices(0, 1) = UBT;
1305
1306 if (debug) {
1307 cerr << "UPSILONB" << endl;
1308 MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
1309 std::string wait;
1310 std::cin >> wait;
1311 }
1312
1313 // CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1314 // CHKERR sub_data->getColIs(&nested_is_cols[0]);
1315 nested_matrices(1, 0) = UPSILONB;
1316
1317 CHKERR DMDestroy(&dm_sub_force);
1318 }
1319
1320 Mat SubA;
1321 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
1322 &sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
1323 &SubA);
1324 nested_matrices(0, 0) = SubA;
1325
1326 CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
1327 CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
1328
1329 if (debug) {
1330 cerr << "Nested SubA" << endl;
1331 MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
1332 }
1333
1334 Mat A;
1335 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
1336 &nested_is_cols[0], &nested_matrices(0, 0), &A);
1337
1338 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1339 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1340
1341 if (debug) {
1342 cerr << "Nested A" << endl;
1343 MatView(A, PETSC_VIEWER_STDOUT_WORLD);
1344 }
1345
1346 Vec D, F;
1347 CHKERR DMCreateGlobalVector(dm_control, &D);
1348 CHKERR DMCreateGlobalVector(dm_control, &F);
1349
1350 // Asseble the right hand side vector
1351 {
1352 CellEngineering::FaceElement face_element(m_field);
1353 face_element.getOpPtrVector().push_back(new OpGetDispX(common_data));
1354 face_element.getOpPtrVector().push_back(new OpGetDispY(common_data));
1355 face_element.getOpPtrVector().push_back(
1356 new OpCell_g(F, eps_u, common_data));
1357 CHKERR DMoFEMLoopFiniteElements(dm_control, "DISPLACEMENTS_PENALTY",
1358 &face_element);
1361 CHKERR VecAssemblyBegin(F);
1362 CHKERR VecAssemblyEnd(F);
1363 }
1364
1365 KSP solver;
1366 // KSP solver;
1367 {
1368 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
1369 CHKERR KSPSetDM(solver, dm_control);
1370 CHKERR KSPSetFromOptions(solver);
1371 CHKERR KSPSetOperators(solver, A, A);
1372 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
1373 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
1374 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
1375 PC pc;
1376 CHKERR KSPGetPC(solver, &pc);
1377 CHKERR PCSetType(pc, PCFIELDSPLIT);
1378 PetscBool is_pcfs = PETSC_FALSE;
1379 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1380 if (is_pcfs) {
1381 CHKERR PCSetOperators(pc, A, A);
1382 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
1383 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
1384 CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1385 CHKERR PCSetUp(pc);
1386 KSP *sub_ksp;
1387 PetscInt n;
1388 CHKERR PCFieldSplitGetSubKSP(pc, &n, &sub_ksp);
1389 {
1390 PC sub_pc_0;
1391 CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
1392 CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
1393 CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
1394 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
1395 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
1396 CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
1397 // CHKERR
1398 // PCFieldSplitSetSchurFactType(sub_pc_0,PC_FIELDSPLIT_SCHUR_FACT_LOWER);
1399 // CHKERR PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_SCHUR);
1400 CHKERR PCSetUp(sub_pc_0);
1401 }
1402 } else {
1403 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1404 "Only works with pre-conditioner PCFIELDSPLIT");
1405 }
1406 CHKERR KSPSetUp(solver);
1407 }
1408
1409 // Solve system of equations
1410 CHKERR KSPSolve(solver, F, D);
1411
1412 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1413 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1414 CHKERR DMoFEMMeshToLocalVector(dm_control, D, INSERT_VALUES,
1415 SCATTER_REVERSE);
1416
1417 if (debug) {
1418 CHKERR VecView(D, PETSC_VIEWER_DRAW_WORLD);
1419 std::string wait;
1420 std::cin >> wait;
1421 }
1422
1423 // Clean sub matrices and sub indices
1424 for (int i = 0; i != 2; i++) {
1425 if (sub_nested_is_rows[i]) {
1426 CHKERR ISDestroy(&sub_nested_is_rows[i]);
1427 }
1428 if (sub_nested_is_cols[i]) {
1429 CHKERR ISDestroy(&sub_nested_is_cols[i]);
1430 }
1431 for (int j = 0; j != 2; j++) {
1432 if (sub_nested_matrices(i, j)) {
1433 CHKERR MatDestroy(&sub_nested_matrices(i, j));
1434 }
1435 }
1436 }
1437 for (int i = 0; i != 2; i++) {
1438 if (nested_is_rows[i]) {
1439 CHKERR ISDestroy(&nested_is_rows[i]);
1440 }
1441 if (nested_is_cols[i]) {
1442 CHKERR ISDestroy(&nested_is_cols[i]);
1443 }
1444 for (int j = 0; j != 2; j++) {
1445 if (nested_matrices(i, j)) {
1446 CHKERR MatDestroy(&nested_matrices(i, j));
1447 }
1448 }
1449 }
1450
1451 CHKERR MatDestroy(&SubA);
1452 CHKERR MatDestroy(&A);
1453 CHKERR VecDestroy(&D);
1454 CHKERR VecDestroy(&F);
1455
1456 CHKERR DMDestroy(&dm_sub_volume_control);
1457
1458 PostProcVolumeOnRefinedMesh post_proc(m_field);
1459 {
1463 // add postpocessing for sresses
1464 post_proc.getOpPtrVector().push_back(new PostProcHookStress(
1465 m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "U",
1466 post_proc.commonData, &elastic.setOfBlocks));
1467 CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC", &post_proc);
1468 CHKERR post_proc.writeFile("out.h5m");
1470 elastic.getLoopFeEnergy().eNergy = 0;
1471 CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC",
1472 &elastic.getLoopFeEnergy());
1473 PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
1474 elastic.getLoopFeEnergy().eNergy);
1475 }
1476
1477 {
1478 PostProcFaceOnRefinedMesh post_proc_face(m_field);
1479 CHKERR post_proc_face.generateReferenceElementMesh();
1480 auto det_ptr = boost::make_shared<VectorDouble>();
1481 auto jac_ptr = boost::make_shared<MatrixDouble>();
1482 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1483 post_proc_face.getOpPtrVector().push_back(new OpCalculateHOJacForFace(jac_ptr));
1484 post_proc_face.getOpPtrVector().push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
1485 if (is_curl) {
1486 post_proc_face.getOpPtrVector().push_back(
1487 new OpSetInvJacHcurlFace(inv_jac_ptr));
1488 post_proc_face.getOpPtrVector().push_back(
1489 new OpVirtualCurlRho("RHO", common_data));
1490 } else {
1491 post_proc_face.getOpPtrVector().push_back(
1492 new OpSetInvJacH1ForFace(inv_jac_ptr));
1493 post_proc_face.getOpPtrVector().push_back(
1494 new OpVirtualPotentialRho("RHO", common_data));
1495 }
1496 post_proc_face.getOpPtrVector().push_back(
1497 new PostProcTraction(m_field, post_proc_face.postProcMesh,
1498 post_proc_face.mapGaussPts, common_data));
1499 CHKERR DMoFEMLoopFiniteElements(dm_control, "D", &post_proc_face);
1500 CHKERR post_proc_face.writeFile("out_tractions.h5m");
1501 }
1502
1503 CHKERR DMDestroy(&dm_control);
1504
1505 }
1507
1509 return 0;
1510}
Implementation of linear elastic material.
int main()
static char help[]
static int debug
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ SIDESET
Definition: definitions.h:147
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
FTensor::Index< 'n', SPACE_DIM > n
@ F
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:219
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:542
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:275
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
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
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
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
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
#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
double D
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
Definition: enable_if.hpp:6
constexpr AssemblyType A
Calculate and assemble g vector.
Definition: CellForces.hpp:268
Calculate and assemble Z matrix.
Definition: CellForces.hpp:393
Calculate and assemble Z matrix.
Definition: CellForces.hpp:316
Calculate and assemble B matrix.
Definition: CellForces.hpp:141
Calculate and assemble D matrix.
Definition: CellForces.hpp:74
Calculate and assemble S matrix.
Definition: CellForces.hpp:202
Post-process tractions.
Definition: CellForces.hpp:529
Shave results on mesh tags for post-processing.
Definition: CellForces.hpp:601
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
Hook equation.
Definition: Hooke.hpp:21
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
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
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
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.
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
Vec & snes_f
residual
Mat & snes_B
preconditioner of jacobian matrix
SNESContext snes_ctx
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
structure grouping operators and data used for calculation of nonlinear elastic element
MyVolumeFE & getLoopFeLhs()
get lhs volume element
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE & getLoopFeRhs()
get rhs volume element
MyVolumeFE & getLoopFeEnergy()
get energy fe
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
Postprocess on face.
MoFEMErrorCode generateReferenceElementMesh()
Operator post-procesing stresses for Hook isotropic material.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
CommonData commonData
Post processing.