v0.14.0
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 
10 In the following code, we solve the problem of identification of forces which
11 result in given displacements on some surface inside the body.
12 
13 In the experiment, a body is covered on top by transparent gel layer, on the
14 interface between movements of markers (bids) is observed. Strain in the body is
15 generated by cells living on the top layer. Here we looking for forces generated
16 by those cells.
17 
18 Following implementing when needed could be easily extended to curved surfaces
19 arbitrary located in the body.
20 
21 \ref mofem_citation
22 
23 \htmlonly
24 <a href="https://doi.org/10.5281/zenodo.439392"><img
25 src="https://zenodo.org/badge/DOI/10.5281/zenodo.439392.svg" alt="DOI"></a> <a
26 href="https://doi.org/10.5281/zenodo.439395"><img
27 src="https://zenodo.org/badge/DOI/10.5281/zenodo.439395.svg" alt="DOI"></a>
28 \endhtmlonly
29 
30 \subsection cell_local Local curl-free formulation
31 
32 This problem can be mathematically described by minimisation problem with the
33 constraints, 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 \\
40 n\cdot\sigma(u) = \rho \quad \textrm{on}\,S_\rho \\
41 \end{array}
42 \right.
43 \f]
44 Applying 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
52 parts \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}
54 +(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,n\cdot\sigma)_{S_\rho}
55 \f]
56 and using second constraint, i.e. static constrain
57 \f[
58 J(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}
61 +(\textrm{grad}[\Upsilon],\sigma)_V - (\Upsilon,\rho)_{S_\rho}.
62 \f]
63 Note that force is enforced in week sense.
64 
65 Notting that
66 \f[
67 \sigma = \mathcal{A} \left( \textrm{grad}[u] \right)^{s}
68 \f]
69 and using symmetry of operator \f$\mathcal{A}\f$,
70 \f[
71 (\textrm{grad}[\Upsilon],\sigma)_V =
72 (\textrm{grad}[\Upsilon],\mathcal{A} \textrm{grad}[u] )_V =
73 (\mathcal{A}^{*}\textrm{grad}[\Upsilon],\textrm{grad}[u] )_V =
74 (\Sigma,\textrm{grad}[u] )_V
75 \f]
76 we finally get
77 \f[
78 J(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}
81 +(\Sigma,\textrm{grad}[u])_V - (\Upsilon,\rho)_{S_\rho}
82 \f]
83 
84 Calculating 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 \\
89 (\sigma,\textrm{grad}[\delta \Upsilon])_V-(\rho,\delta \Upsilon)_{S_\rho} &= 0
90 \\
91 (\epsilon_\rho \rho,\delta \rho)_{S_\rho}-(\Upsilon,\delta \rho)_{S_\rho} &= 0
92 \end{array}
93 \right.
94 \f]
95 and applying finite element discretization, above equations are expressed
96 by linear algebraic equations suitable for computer calculations
97 \f[
98 \left[
99 \begin{array}{ccc}
100 S & K_{u \Upsilon} & 0 \\
101 K_{\Upsilon u} & 0 & -B^\textrm{T} \\
102 0 & -B & D
103 \end{array}
104 \right]
105 \left[
106 \begin{array}{c}
107 u \\ \Upsilon \\ \rho
108 \end{array}
109 \right]
110 =
111 \left[
112 \begin{array}{c}
113 S u_d \\ 0 \\ 0
114 \end{array}
115 \right]
116 \f]
117 and finally swapping first two rows, we finally get
118 \f[
119 \left[
120 \begin{array}{ccc}
121 K & 0 & -B^\textrm{T} \\
122 S & K & 0 \\
123 0 & -B & D
124 \end{array}
125 \right]
126 \left[
127 \begin{array}{c}
128 u \\ \Upsilon \\ \rho
129 \end{array}
130 \right]
131 =
132 \left[
133 \begin{array}{c}
134 0 \\ S u_d \\ 0
135 \end{array}
136 \right]
137 \f]
138 where
139 \f[
140 S=\epsilon_u^{-1} I\\
141 \epsilon_\rho^\textrm{abs} = \epsilon_\rho^\text{ref}\epsilon_u^{-1}
142 \f]
143 
144 Displacement field is \f$\mathbf{u}\f$ and \f$\Upsilon\f$, force field on top
145 surface is \f$\boldsymbol\rho\f$.
146 
147 We not yet explored nature of cell forces. Since the cells are attached to the
148 body surface and anything else, the body as results those forces can not be
149 subjected to rigid body motion.
150 
151 We assume that forces are generated by 2d objects, as results only tangential
152 forces can be produced by those forces if the surface is planar. For non-planar
153 surfaces addition force normal to the surface is present, although the magnitude
154 of this force could be calculated purely from geometrical considerations, thus
155 additional physical equations are not needed.
156 
157 Assuming that straight and very short pre-stressed fibres generate cell forces;
158 a force field is curl-free. Utilising that forces can be expressed by potential
159 field as follows
160 \f[
161 \rho = \frac{\partial \Phi_\rho}{\partial \mathbf{x}}
162 \f]
163 where \f$\Phi\f$ is force potential field.
164 
165 \subsubsection cell_local_approx Approximation
166 
167 For 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
169 potential field \f$\Phi\f$ which derivatives gives potential forces.
170 
171 \subsection cell_nonlocal Nonlocal, i.e. weak curl-free
172 
173 In this variant generalisation of the local model to a weakly enforced force
174 curl-free cell force field is considered. This generalisation is a consequence
175 of the observation that cell has some small but finite size. Moreover, a cell
176 has a complex pre-stressed structure which can transfer shear forces.
177 
178 We 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 \\
187 n\cdot\sigma = \rho
188 \end{array}
189 \right.
190 \f]
191 and 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 \\
200 n\cdot\sigma = \rho
201 \end{array}
202 \right.
203 \f]
204 where parameter \f$\epsilon_l\f$ controls curl-free term, if this parameter
205 approach zero, this formulation converge to local variant presented above.
206 
207 Applying 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
216 give following system of equations \f[ \left[ \begin{array}{ccc}
217 S & K_{u \Upsilon} & 0 \\
218 K_{\Upsilon u} & 0 & -B^\textrm{T} \\
219 0 & -B & D
220 \end{array}
221 \right]
222 \left[
223 \begin{array}{c}
224 u \\ \Upsilon \\ \rho
225 \end{array}
226 \right]
227 =
228 \left[
229 \begin{array}{c}
230 S u_d \\ 0 \\ 0
231 \end{array}
232 \right]
233 \f]
234 and swapping the first two rows we get
235 \f[
236 \left[
237 \begin{array}{ccc}
238 K & 0 & -B^\textrm{T} \\
239 S & K & 0 \\
240 0 & -B & D
241 \end{array}
242 \right]
243 \left[
244 \begin{array}{c}
245 u \\ \Upsilon \\ \rho
246 \end{array}
247 \right]
248 =
249 \left[
250 \begin{array}{c}
251 0 \\ S u_d \\ 0
252 \end{array}
253 \right]
254 \f]
255 
256 \subsubsection cell_nonlocal_approx Approximation
257 
258 For 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 
264 In this implementation, the block structure of the matrix is utilised. The \e
265 PCFIELDSPLIT
266 (see
267 <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html>)
268 is 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 
272 The sub-netsetd matrix is created containing two upper-left blocks,
273 \f[
274 X = \left[
275 \begin{array}{cc}
276 K & 0 \\
277 S & K
278 \end{array}
279 \right]
280 \f]
281 this netsed matrix is blocked in
282 \f[
283 A= \left[
284 \begin{array}{cc}
285 X & V \\
286 H & D
287 \end{array}
288 \right]
289 \f]
290 where
291 \f[
292 V= \left[
293 \begin{array}{c}
294 -B^\textrm{T} \\
295 0
296 \end{array}
297 \right]
298 \f]
299 and
300 \f[
301 H= \left[
302 \begin{array}{cc}
303 0 & -B^\textrm{T}
304 \end{array}
305 \right]
306 \f]
307 
308 The system associated with X matrix is solved using \e PCFIELDSPLIT and
309 multiplicative relaxation scheme. In following implementation sub-matrix, K is
310 factorised only once. The system related to matrix \f$A\f$ composed form
311 matrices \f$X,D,V,H\f$ is solved using Schur complement. Note that here we
312 exploit MoFEM capability to create sub-problems in easy way.
313 
314 \section cell_running_code Running code
315 
316 Approximations 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 
326 If only one layer is specified, another thin polymer layer could be created
327 automatically 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 
338 Config file
339 \include users_modules/cell_engineering/examples/block_config.in
340 
341 Calculating forces
342 \code
343 mpirun -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 
364 As 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
376 metric tensors and taking into account curvature.
377 
378 \todo Instead of elastic material we can use GEL model developed in other
379 module. That will allow to take into account drying and other rheological
380 effects.
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
385 functional 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
392  * the terms of the GNU Lesser General Public License as published by the
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
402  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
403 
404 #include <BasicFiniteElements.hpp>
405 using namespace MoFEM;
406 #include <Hooke.hpp>
407 #include <CellForces.hpp>
408 
409 static char help[] = "-my_block_config set block data\n"
410  "\n";
411 
412 namespace CellEngineering {
413 
415  int oRder;
416  double yOung;
417  double pOisson;
418  BlockOptionData() : oRder(-1), yOung(-1), pOisson(-2) {}
419 };
420 
421 } // namespace CellEngineering
422 
423 using namespace boost::numeric;
424 using namespace CellEngineering;
425 
426 #include <boost/program_options.hpp>
427 using namespace std;
428 namespace po = boost::program_options;
429 
430 static int debug = 0;
431 
432 int 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;
487  option = "PARALLEL=READ_PART;"
488  "PARALLEL_RESOLVE_SHARED_ENTS;"
489  "PARTITION=PARALLEL_PARTITION;";
490  // option = "";
491  CHKERR moab.load_file(mesh_file_name, 0, option);
492 
493  // Create MoFEM (Joseph) database
494  MoFEM::Core core(moab);
495  MoFEM::Interface &m_field = core;
496 
497  MeshsetsManager *mmanager_ptr;
498  CHKERR m_field.query_interface(mmanager_ptr);
499  // print bcs
500  CHKERR mmanager_ptr->printDisplacementSet();
501  CHKERR mmanager_ptr->printForceSet();
502  // print block sets with materials
503  CHKERR mmanager_ptr->printMaterialsSet();
504 
505  // stl::bitset see for more details
506  BitRefLevel bit_level0;
507  bit_level0.set(0);
508  {
509  Range ents3d;
510  CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
511  CHKERR m_field.seed_ref_level(ents3d, bit_level0, false);
512  }
513 
514  // set app. order
515 
516  std::vector<Range> setOrderToEnts(10);
517 
518  // configure blocks by parsing config file
519  // it allow to set approximation order for each block independently
520  Range set_order_ents;
521  std::map<int, BlockOptionData> block_data;
522  if (flg_block_config) {
523  double read_eps_u, read_eps_rho, read_eps_l;
524  try {
525  ifstream ini_file(block_config_file);
526  if (!ini_file.is_open()) {
527  SETERRQ(PETSC_COMM_SELF, 1,
528  "*** -my_block_config does not exist ***");
529  }
530  // std::cerr << block_config_file << std::endl;
531  po::variables_map vm;
532  po::options_description config_file_options;
533  config_file_options.add_options()(
534  "eps_u", po::value<double>(&read_eps_u)->default_value(-1))(
535  "eps_rho", po::value<double>(&read_eps_rho)->default_value(-1))(
536  "eps_l", po::value<double>(&read_eps_l)->default_value(-1));
538  std::ostringstream str_order;
539  str_order << "block_" << it->getMeshsetId() << ".displacement_order";
540  config_file_options.add_options()(
541  str_order.str().c_str(),
542  po::value<int>(&block_data[it->getMeshsetId()].oRder)
543  ->default_value(order));
544  std::ostringstream str_cond;
545  str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
546  config_file_options.add_options()(
547  str_cond.str().c_str(),
548  po::value<double>(&block_data[it->getMeshsetId()].yOung)
549  ->default_value(-1));
550  std::ostringstream str_capa;
551  str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
552  config_file_options.add_options()(
553  str_capa.str().c_str(),
554  po::value<double>(&block_data[it->getMeshsetId()].pOisson)
555  ->default_value(-2));
556  }
557  po::parsed_options parsed =
558  parse_config_file(ini_file, config_file_options, true);
559  store(parsed, vm);
560  po::notify(vm);
562  if (block_data[it->getMeshsetId()].oRder == -1)
563  continue;
564  if (block_data[it->getMeshsetId()].oRder == order)
565  continue;
566  PetscPrintf(PETSC_COMM_WORLD, "Set block %d order to %d\n",
567  it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
568  Range block_ents;
569  CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
570  true);
571  // block_ents = block_ents.subset_by_type(MBTET);
572  Range nodes;
573  CHKERR moab.get_connectivity(block_ents, nodes, true);
574  Range ents_to_set_order, ents3d;
575  CHKERR moab.get_adjacencies(nodes, 3, false, ents3d,
576  moab::Interface::UNION);
577  CHKERR moab.get_adjacencies(ents3d, 2, false, ents_to_set_order,
578  moab::Interface::UNION);
579  CHKERR moab.get_adjacencies(ents3d, 1, false, ents_to_set_order,
580  moab::Interface::UNION);
581  ents_to_set_order = subtract(
582  ents_to_set_order, ents_to_set_order.subset_by_type(MBQUAD));
583  ents_to_set_order = subtract(
584  ents_to_set_order, ents_to_set_order.subset_by_type(MBPRISM));
585  set_order_ents.merge(ents3d);
586  set_order_ents.merge(ents_to_set_order);
587  setOrderToEnts[block_data[it->getMeshsetId()].oRder].merge(
588  set_order_ents);
589  }
590  CHKERR m_field.synchronise_entities(set_order_ents, 0);
591  std::vector<std::string> additional_parameters;
592  additional_parameters =
593  collect_unrecognized(parsed.options, po::include_positional);
594  for (std::vector<std::string>::iterator vit =
595  additional_parameters.begin();
596  vit != additional_parameters.end(); vit++) {
597  CHKERR PetscPrintf(PETSC_COMM_WORLD,
598  "** WARNING Unrecognized option %s\n",
599  vit->c_str());
600  }
601  } catch (const std::exception &ex) {
602  std::ostringstream ss;
603  ss << ex.what() << std::endl;
604  SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
605  }
606  if (read_eps_u > 0) {
607  eps_u = read_eps_u;
608  };
609  if (read_eps_rho > 0) {
610  eps_rho = read_eps_rho;
611  }
612  if (read_eps_l > 0) {
613  eps_l = read_eps_l;
614  }
615  }
616 
617  PetscPrintf(PETSC_COMM_WORLD, "epsU = %6.4e epsRho = %6.4e\n", eps_u,
618  eps_rho);
619 
620  // Fields
621  CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 3, MB_TAG_SPARSE,
622  MF_ZERO);
623  CHKERR m_field.add_field("UPSILON", H1, AINSWORTH_LEGENDRE_BASE, 3,
624  MB_TAG_SPARSE, MF_ZERO);
625  if (is_curl) {
626  CHKERR m_field.add_field("RHO", HCURL, AINSWORTH_LEGENDRE_BASE, 1);
627  } else {
628  CHKERR m_field.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1);
629  }
630 
631  // add entitities (by tets) to the field
632  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "U");
633  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "UPSILON");
634  CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "U");
635  CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "UPSILON");
636 
637  CHKERR m_field.synchronise_field_entities("U");
638  CHKERR m_field.synchronise_field_entities("UPSILON");
639  CHKERR m_field.synchronise_field_entities("RHO");
640 
641  Range vertex_to_fix;
642  Range edges_to_fix;
643  Range ents_1st_layer;
644  // If problem is partitioned meshset culd not exist on this part
645  if (mmanager_ptr->checkMeshset(202, SIDESET)) {
646  CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 2,
647  ents_1st_layer, true);
648  CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 0,
649  vertex_to_fix, false);
650  CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 1, edges_to_fix,
651  false);
652  if (vertex_to_fix.size() != 1 && !vertex_to_fix.empty()) {
653  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
654  "Should be one vertex only, but is %d", vertex_to_fix.size());
655  }
656  }
657  CHKERR m_field.synchronise_entities(ents_1st_layer, 0);
659  ents_1st_layer.subset_by_type(MBTRI), MBTRI, "RHO");
660  Range ents_2nd_layer;
661  // If problem is partitioned meshset culd not exist on this part
662  if (mmanager_ptr->checkMeshset(101, SIDESET)) {
663  CHKERR mmanager_ptr->getEntitiesByDimension(101, SIDESET, 2,
664  ents_2nd_layer, true);
665  }
666  CHKERR m_field.synchronise_entities(ents_2nd_layer, 0);
667 
668  for (int oo = 2; oo != setOrderToEnts.size(); oo++) {
669  if (setOrderToEnts[oo].size() > 0) {
670  CHKERR m_field.synchronise_entities(setOrderToEnts[oo], 0);
671  CHKERR m_field.set_field_order(setOrderToEnts[oo], "U", oo);
672  CHKERR m_field.set_field_order(setOrderToEnts[oo], "UPSILON", oo);
673  }
674  }
675 
676  const int through_thickness_order = 2;
677  {
678  Range ents3d;
679  CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
680  Range ents;
681  CHKERR moab.get_adjacencies(ents3d, 2, false, ents,
682  moab::Interface::UNION);
683  CHKERR moab.get_adjacencies(ents3d, 1, false, ents,
684  moab::Interface::UNION);
685 
686  Range prisms;
687  CHKERR moab.get_entities_by_type(0, MBPRISM, prisms);
688  {
689  Range quads;
690  CHKERR moab.get_adjacencies(prisms, 2, false, quads,
691  moab::Interface::UNION);
692  Range prism_tris;
693  prism_tris = quads.subset_by_type(MBTRI);
694  quads = subtract(quads, prism_tris);
695  Range quads_edges;
696  CHKERR moab.get_adjacencies(quads, 1, false, quads_edges,
697  moab::Interface::UNION);
698  Range prism_tris_edges;
699  CHKERR moab.get_adjacencies(prism_tris, 1, false, prism_tris_edges,
700  moab::Interface::UNION);
701  quads_edges = subtract(quads_edges, prism_tris_edges);
702  prisms.merge(quads);
703  prisms.merge(quads_edges);
704  }
705 
706  ents.merge(ents3d);
707  ents = subtract(ents, set_order_ents);
708  ents = subtract(ents, prisms);
709 
710  CHKERR m_field.synchronise_entities(ents, 0);
711  CHKERR m_field.synchronise_entities(prisms, 0);
712 
713  CHKERR m_field.set_field_order(ents, "U", order);
714  CHKERR m_field.set_field_order(ents, "UPSILON", order);
715  // approx. order through thickness to 2
716  CHKERR m_field.set_field_order(prisms, "U", through_thickness_order);
717  CHKERR m_field.set_field_order(prisms, "UPSILON",
718  through_thickness_order);
719  }
720  CHKERR m_field.set_field_order(0, MBVERTEX, "U", 1);
721  CHKERR m_field.set_field_order(0, MBVERTEX, "UPSILON", 1);
722 
723  if (is_curl) {
724  CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
725  CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
726  } else {
727  CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
728  CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
729  CHKERR m_field.set_field_order(0, MBVERTEX, "RHO", 1);
730  }
731 
732  // Add elastic element
733  boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(new Hooke<adouble>());
734  boost::shared_ptr<Hooke<double>> hooke_double_ptr(new Hooke<double>());
735  NonlinearElasticElement elastic(m_field, 2);
736  CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
737  CHKERR elastic.addElement("ELASTIC", "U");
738  CHKERR elastic.setOperators("U", "MESH_NODE_POSITIONS", false, true);
739  // Add prisms
740  CellEngineering::FatPrism fat_prism_rhs(m_field);
741  CellEngineering::FatPrism fat_prism_lhs(m_field);
742  {
743  CHKERR m_field.add_finite_element("ELASTIC_PRISM", MF_ZERO);
744  CHKERR m_field.modify_finite_element_add_field_row("ELASTIC_PRISM", "U");
745  CHKERR m_field.modify_finite_element_add_field_col("ELASTIC_PRISM", "U");
746  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC_PRISM", "U");
747  Range ents_2nd_layer;
748  CHKERR mmanager_ptr->getEntitiesByDimension(2, BLOCKSET, 3,
749  elastic.setOfBlocks[2].tEts);
751  elastic.setOfBlocks[2].tEts, MBPRISM, "ELASTIC_PRISM");
752  MatrixDouble inv_jac;
753  // right hand side operators
754  fat_prism_rhs.getOpPtrVector().push_back(
756  fat_prism_rhs.getOpPtrVector().push_back(
757  new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac));
758  fat_prism_rhs.getOpPtrVector().push_back(
760  "U", elastic.commonData));
761  fat_prism_rhs.getOpPtrVector().push_back(
763  "U", elastic.setOfBlocks[2], elastic.commonData, 2, false, false,
764  true));
765  fat_prism_rhs.getOpPtrVector().push_back(
767  "U", elastic.setOfBlocks[2], elastic.commonData));
768  // Left hand side operators
769  fat_prism_lhs.getOpPtrVector().push_back(
771  fat_prism_lhs.getOpPtrVector().push_back(
772  new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac));
773  fat_prism_lhs.getOpPtrVector().push_back(
775  "U", elastic.commonData));
776  fat_prism_lhs.getOpPtrVector().push_back(
778  "U", elastic.setOfBlocks[2], elastic.commonData, 2, true, false,
779  true));
780  fat_prism_lhs.getOpPtrVector().push_back(
782  "U", "U", elastic.setOfBlocks[2], elastic.commonData));
783  }
784 
785  // Update material parameters
787  int id = it->getMeshsetId();
788  if (block_data[id].yOung > 0) {
789  elastic.setOfBlocks[id].E = block_data[id].yOung;
790  CHKERR PetscPrintf(PETSC_COMM_WORLD,
791  "Block %d set Young modulus %3.4g\n", id,
792  elastic.setOfBlocks[id].E);
793  }
794  if (block_data[id].pOisson >= -1) {
795  elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
796  CHKERR PetscPrintf(PETSC_COMM_WORLD,
797  "Block %d set Poisson ratio %3.4g\n", id,
798  elastic.setOfBlocks[id].PoissonRatio);
799  }
800  }
801 
802  // build field
803  CHKERR m_field.build_fields();
804 
805  // Control elements
806  CHKERR m_field.add_finite_element("KUPSUPS");
807  CHKERR m_field.modify_finite_element_add_field_row("KUPSUPS", "UPSILON");
808  CHKERR m_field.modify_finite_element_add_field_col("KUPSUPS", "UPSILON");
809  CHKERR m_field.modify_finite_element_add_field_data("KUPSUPS", "UPSILON");
810  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTET, "KUPSUPS");
811  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBPRISM, "KUPSUPS");
812 
813  CHKERR m_field.add_finite_element("DISPLACEMENTS_PENALTY");
814  CHKERR m_field.modify_finite_element_add_field_row("DISPLACEMENTS_PENALTY",
815  "UPSILON");
816  CHKERR m_field.modify_finite_element_add_field_col("DISPLACEMENTS_PENALTY",
817  "U");
818  CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
819  "UPSILON");
820  CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
821  "U");
822  CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
823  "DISP_X");
824  CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
825  "DISP_Y");
826  CHKERR m_field.add_ents_to_finite_element_by_type(ents_2nd_layer, MBTRI,
827  "DISPLACEMENTS_PENALTY");
828 
829  // Add element to calculate residual on 1st layer
830  CHKERR m_field.add_finite_element("BT");
831  CHKERR m_field.modify_finite_element_add_field_row("BT", "U");
832  CHKERR m_field.modify_finite_element_add_field_row("BT", "UPSILON");
833  CHKERR m_field.modify_finite_element_add_field_col("BT", "RHO");
834  CHKERR m_field.modify_finite_element_add_field_data("BT", "U");
835  CHKERR m_field.modify_finite_element_add_field_data("BT", "UPSILON");
836  CHKERR m_field.modify_finite_element_add_field_data("BT", "RHO");
837  CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
838  "BT");
839 
840  CHKERR m_field.add_finite_element("B");
841  CHKERR m_field.modify_finite_element_add_field_row("B", "RHO");
843  CHKERR m_field.modify_finite_element_add_field_col("B", "UPSILON");
845  CHKERR m_field.modify_finite_element_add_field_data("B", "UPSILON");
846  CHKERR m_field.modify_finite_element_add_field_data("B", "RHO");
847  CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
848  "B");
849 
850  // Add element to calculate residual on 1st layer
851  CHKERR m_field.add_finite_element("D");
852  CHKERR m_field.modify_finite_element_add_field_row("D", "RHO");
853  CHKERR m_field.modify_finite_element_add_field_col("D", "RHO");
854  CHKERR m_field.modify_finite_element_add_field_data("D", "RHO");
855  CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
856  "D");
857 
858  // build finite elemnts
859  CHKERR m_field.build_finite_elements();
860  // build adjacencies
861  CHKERR m_field.build_adjacencies(bit_level0);
862 
863  // Register MOFEM DM
864  DMType dm_name = "MOFEM";
865  CHKERR DMRegister_MoFEM(dm_name);
866 
867  DM dm_control;
868  {
869  // Craete dm_control instance
870  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
871  CHKERR DMSetType(dm_control, dm_name);
872  // set dm_control data structure which created mofem datastructures
873  CHKERR DMMoFEMCreateMoFEM(dm_control, &m_field, "CONTROL_PROB",
874  bit_level0);
875  CHKERR DMSetFromOptions(dm_control);
876  CHKERR DMMoFEMSetSquareProblem(dm_control, PETSC_TRUE);
877  CHKERR DMMoFEMSetIsPartitioned(dm_control, PETSC_TRUE);
878  // add elements to dm_control
879  CHKERR DMMoFEMAddElement(dm_control, "ELASTIC");
880  CHKERR DMMoFEMAddElement(dm_control, "ELASTIC_PRISM");
881  CHKERR DMMoFEMAddElement(dm_control, "KUPSUPS");
882  CHKERR DMMoFEMAddElement(dm_control, "DISPLACEMENTS_PENALTY");
883  CHKERR DMMoFEMAddElement(dm_control, "B");
884  CHKERR DMMoFEMAddElement(dm_control, "BT");
885  CHKERR DMMoFEMAddElement(dm_control, "D");
886  CHKERR DMSetUp(dm_control);
887  }
888 
889  MatrixDouble inv_jac;
890  CellEngineering::CommonData common_data;
891 
892  ublas::matrix<Mat> nested_matrices(2, 2);
893  ublas::vector<IS> nested_is_rows(2);
894  ublas::vector<IS> nested_is_cols(2);
895  for (int i = 0; i != 2; i++) {
896  nested_is_rows[i] = PETSC_NULL;
897  nested_is_cols[i] = PETSC_NULL;
898  for (int j = 0; j != 2; j++) {
899  nested_matrices(i, j) = PETSC_NULL;
900  }
901  }
902 
903  ublas::matrix<Mat> sub_nested_matrices(2, 2);
904  ublas::vector<IS> sub_nested_is_rows(2);
905  ublas::vector<IS> sub_nested_is_cols(2);
906  for (int i = 0; i != 2; i++) {
907  sub_nested_is_rows[i] = PETSC_NULL;
908  sub_nested_is_cols[i] = PETSC_NULL;
909  for (int j = 0; j != 2; j++) {
910  sub_nested_matrices(i, j) = PETSC_NULL;
911  }
912  }
913 
914  DM dm_sub_volume_control;
915  {
916  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_volume_control);
917  CHKERR DMSetType(dm_sub_volume_control, dm_name);
918  // set dm_sub_volume_control data structure which created mofem data
919  // structures
920  CHKERR DMMoFEMCreateSubDM(dm_sub_volume_control, dm_control,
921  "SUB_CONTROL_PROB");
922  CHKERR DMMoFEMSetSquareProblem(dm_sub_volume_control, PETSC_TRUE);
923  CHKERR DMMoFEMAddSubFieldRow(dm_sub_volume_control, "U");
924  CHKERR DMMoFEMAddSubFieldRow(dm_sub_volume_control, "UPSILON");
925  CHKERR DMMoFEMAddSubFieldCol(dm_sub_volume_control, "U");
926  CHKERR DMMoFEMAddSubFieldCol(dm_sub_volume_control, "UPSILON");
927  // add elements to dm_sub_volume_control
928  CHKERR DMSetUp(dm_sub_volume_control);
929 
930  const Problem *prb_ptr;
931  CHKERR m_field.get_problem("SUB_CONTROL_PROB", &prb_ptr);
932  boost::shared_ptr<Problem::SubProblemData> sub_data =
933  prb_ptr->getSubData();
934 
935  CHKERR sub_data->getRowIs(&nested_is_rows[0]);
936  CHKERR sub_data->getColIs(&nested_is_cols[0]);
937  // That will be filled at the end
938  nested_matrices(0, 0) = PETSC_NULL;
939  }
940 
941  {
942  DM dm_sub_sub_elastic;
943 
944  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_sub_elastic);
945  CHKERR DMSetType(dm_sub_sub_elastic, dm_name);
946  // set dm_sub_sub_elastic data structure which created mofem data
947  // structures
948  CHKERR DMMoFEMCreateSubDM(dm_sub_sub_elastic, dm_sub_volume_control,
949  "ELASTIC_PROB");
950  CHKERR DMMoFEMSetSquareProblem(dm_sub_sub_elastic, PETSC_TRUE);
951  CHKERR DMMoFEMAddElement(dm_sub_sub_elastic, "ELASTIC");
952  CHKERR DMMoFEMAddElement(dm_sub_sub_elastic, "ELASTIC_PRISM");
953  CHKERR DMMoFEMAddSubFieldRow(dm_sub_sub_elastic, "U");
954  CHKERR DMMoFEMAddSubFieldCol(dm_sub_sub_elastic, "U");
955  // add elements to dm_sub_sub_elastic
956  CHKERR DMSetUp(dm_sub_sub_elastic);
957 
958  Mat Kuu;
959  Vec Du, Fu;
960  CHKERR DMCreateMatrix(dm_sub_sub_elastic, &Kuu);
961  CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Du);
962  CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Fu);
963  CHKERR MatZeroEntries(Kuu);
964  CHKERR VecZeroEntries(Du);
965  CHKERR VecZeroEntries(Fu);
966  CHKERR DMoFEMMeshToLocalVector(dm_sub_sub_elastic, Du, INSERT_VALUES,
967  SCATTER_REVERSE);
968 
969  // Manage Dirichlet bC
970  boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
971  dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
972  new DirichletDisplacementBc(m_field, "U", Kuu, Du, Fu));
973  dirihlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
974  dirihlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
975  // preproc
976  CHKERR DMoFEMPreProcessFiniteElements(dm_sub_sub_elastic,
977  dirihlet_bc_ptr.get());
978  // internal force vector (to take into account Dirichlet boundary
979  // conditions)
980  elastic.getLoopFeRhs().snes_f = Fu;
981  fat_prism_rhs.snes_f = Fu;
982  CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
983  &elastic.getLoopFeRhs());
984  CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
985  &fat_prism_rhs);
986  // elastic element matrix
987  elastic.getLoopFeLhs().snes_B = Kuu;
988  fat_prism_lhs.snes_B = Kuu;
989  CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
990  &elastic.getLoopFeLhs());
991  CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
992  &fat_prism_lhs);
993  // postproc
994  CHKERR DMoFEMPostProcessFiniteElements(dm_sub_sub_elastic,
995  dirihlet_bc_ptr.get());
996  CHKERR VecGhostUpdateBegin(Fu, ADD_VALUES, SCATTER_REVERSE);
997  CHKERR VecGhostUpdateEnd(Fu, ADD_VALUES, SCATTER_REVERSE);
998  CHKERR VecAssemblyBegin(Fu);
999  CHKERR VecAssemblyEnd(Fu);
1000  CHKERR VecScale(Fu, -1);
1001  CHKERR VecDestroy(&Du);
1002  CHKERR VecDestroy(&Fu);
1003 
1004  const Problem *prb_ptr;
1005  CHKERR m_field.get_problem("ELASTIC_PROB", &prb_ptr);
1006  boost::shared_ptr<Problem::SubProblemData> sub_data =
1007  prb_ptr->getSubData();
1008 
1009  CHKERR sub_data->getRowIs(&sub_nested_is_rows[0]);
1010  CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1011  sub_nested_matrices(0, 0) = Kuu;
1012  IS isUpsilon;
1013  CHKERR m_field.query_interface<ISManager>()
1014  ->isCreateFromProblemFieldToOtherProblemField(
1015  "ELASTIC_PROB", "U", ROW, "SUB_CONTROL_PROB", "UPSILON", ROW,
1016  PETSC_NULL, &isUpsilon);
1017  sub_nested_is_rows[1] = isUpsilon;
1018  sub_nested_is_cols[1] = isUpsilon;
1019  sub_nested_matrices(1, 1) = Kuu;
1020  PetscObjectReference((PetscObject)Kuu);
1021  PetscObjectReference((PetscObject)isUpsilon);
1022 
1023  // Matrix View
1024  if (debug) {
1025  cerr << "Kuu" << endl;
1026  MatView(Kuu, PETSC_VIEWER_DRAW_WORLD);
1027  std::string wait;
1028  std::cin >> wait;
1029  }
1030 
1031  CHKERR DMDestroy(&dm_sub_sub_elastic);
1032  }
1033 
1034  {
1035  DM dm_sub_disp_penalty;
1036 
1037  // Craete dm_control instance
1038  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_disp_penalty);
1039  CHKERR DMSetType(dm_sub_disp_penalty, dm_name);
1040  // set dm_sub_disp_penalty data structure which created mofem
1041  // datastructures
1042  CHKERR DMMoFEMCreateSubDM(dm_sub_disp_penalty, dm_sub_volume_control,
1043  "S_PROB");
1044  CHKERR DMMoFEMSetSquareProblem(dm_sub_disp_penalty, PETSC_FALSE);
1045  CHKERR DMMoFEMAddSubFieldRow(dm_sub_disp_penalty, "UPSILON");
1046  CHKERR DMMoFEMAddSubFieldCol(dm_sub_disp_penalty, "U");
1047  // add elements to dm_sub_disp_penalty
1048  CHKERR DMMoFEMAddElement(dm_sub_disp_penalty, "DISPLACEMENTS_PENALTY");
1049  CHKERR DMSetUp(dm_sub_disp_penalty);
1050 
1051  Mat S;
1052  CHKERR DMCreateMatrix(dm_sub_disp_penalty, &S);
1053  CHKERR MatZeroEntries(S);
1054 
1055  CellEngineering::FaceElement face_element(m_field);
1056  face_element.getOpPtrVector().push_back(new OpCellS(S, eps_u));
1057  CHKERR DMoFEMLoopFiniteElements(dm_sub_disp_penalty,
1058  "DISPLACEMENTS_PENALTY", &face_element);
1059  CHKERR MatAssemblyBegin(S, MAT_FLUSH_ASSEMBLY);
1060  CHKERR MatAssemblyEnd(S, MAT_FLUSH_ASSEMBLY);
1061 
1062  // // Matrix View
1063  if (debug) {
1064  cerr << "S" << endl;
1065  MatView(S, PETSC_VIEWER_DRAW_WORLD);
1066  std::string wait;
1067  std::cin >> wait;
1068  }
1069 
1070  const Problem *problem_ptr;
1071  CHKERR m_field.get_problem("S_PROB", &problem_ptr);
1072  boost::shared_ptr<Problem::SubProblemData> sub_data =
1073  problem_ptr->getSubData();
1074  // CHKERR sub_data->getRowIs(&sub_nested_is_rows[1]);
1075  // CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1076  sub_nested_matrices(1, 0) = S;
1077 
1078  CHKERR DMDestroy(&dm_sub_disp_penalty);
1079  }
1080 
1081  // Calculate penalty matrix
1082  {
1083  DM dm_sub_force_penalty;
1084 
1085  // Craete dm_control instance
1086  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force_penalty);
1087  CHKERR DMSetType(dm_sub_force_penalty, dm_name);
1088  // set dm_sub_force_penalty data structure which created mofem
1089  // datastructures
1090  CHKERR DMMoFEMCreateSubDM(dm_sub_force_penalty, dm_control, "D_PROB");
1091  CHKERR DMMoFEMSetSquareProblem(dm_sub_force_penalty, PETSC_TRUE);
1092  CHKERR DMMoFEMAddSubFieldRow(dm_sub_force_penalty, "RHO");
1093  CHKERR DMMoFEMAddSubFieldCol(dm_sub_force_penalty, "RHO");
1094  // add elements to dm_sub_force_penalty
1095  CHKERR DMMoFEMAddElement(dm_sub_force_penalty, "D");
1096  CHKERR DMSetUp(dm_sub_force_penalty);
1097 
1098  Mat D;
1099  CHKERR DMCreateMatrix(dm_sub_force_penalty, &D);
1100  CHKERR MatZeroEntries(D);
1101 
1102  {
1103  CellEngineering::FaceElement face_d_matrix(m_field);
1104  face_d_matrix.getOpPtrVector().push_back(
1105  new OpCalculateInvJacForFace(inv_jac));
1106  if (is_curl) {
1107  face_d_matrix.getOpPtrVector().push_back(
1108  new OpSetInvJacHcurlFace(inv_jac));
1109  face_d_matrix.getOpPtrVector().push_back(
1110  new OpCellCurlD(D, eps_rho / eps_u, eps_l * eps_u));
1111  } else {
1112  face_d_matrix.getOpPtrVector().push_back(
1113  new OpSetInvJacH1ForFace(inv_jac));
1114  face_d_matrix.getOpPtrVector().push_back(
1115  new OpCellPotentialD(D, eps_rho / eps_u));
1116  }
1117  CHKERR DMoFEMLoopFiniteElements(dm_sub_force_penalty, "D",
1118  &face_d_matrix);
1119  }
1120  CHKERR MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY);
1121  CHKERR MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY);
1122 
1123  const Problem *problem_ptr;
1124  CHKERR m_field.get_problem("D_PROB", &problem_ptr);
1125 
1126  // Zero rows, force field is given by gradients of potential field, so one
1127  // of the values has to be fixed like for rigid body motion.
1128  if (is_curl == PETSC_FALSE) {
1129  int nb_dofs_to_fix = 0;
1130  int index_to_fix = 0;
1131  if (!vertex_to_fix.empty()) {
1132  boost::shared_ptr<NumeredDofEntity> dof_ptr;
1133  CHKERR problem_ptr->getDofByNameEntAndEntDofIdx(
1134  m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
1135  dof_ptr);
1136  if (dof_ptr) {
1137  if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1138  nb_dofs_to_fix = 1;
1139  index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1140  cerr << *dof_ptr << endl;
1141  }
1142  }
1143  }
1144  CHKERR MatZeroRowsColumns(D, nb_dofs_to_fix, &index_to_fix,
1145  eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1146  } else {
1147  std::vector<int> dofs_to_fix;
1148  for (auto p_eit = edges_to_fix.pair_begin();
1149  p_eit != edges_to_fix.pair_end(); ++p_eit) {
1150  auto bit_number = m_field.get_field_bit_number("RHO");
1151  auto row_dofs = problem_ptr->numeredRowDofsPtr;
1152  auto lo = row_dofs->lower_bound(
1153  FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
1154  auto hi =
1155  row_dofs->lower_bound(FieldEntity::getHiLocalEntityBitNumber(
1156  bit_number, p_eit->second));
1157  for (; lo != hi; ++lo)
1158  dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1159  }
1160  CHKERR MatZeroRowsColumns(D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1161  eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1162  }
1163 
1164  // Matrix View
1165  if (debug) {
1166  cerr << "D" << endl;
1167  MatView(D, PETSC_VIEWER_DRAW_WORLD);
1168  std::string wait;
1169  std::cin >> wait;
1170  }
1171 
1172  boost::shared_ptr<Problem::SubProblemData> sub_data =
1173  problem_ptr->getSubData();
1174  CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1175  CHKERR sub_data->getColIs(&nested_is_cols[1]);
1176  nested_matrices(1, 1) = D;
1177 
1178  CHKERR DMDestroy(&dm_sub_force_penalty);
1179  }
1180 
1181  {
1182  DM dm_sub_force;
1183 
1184  // Craete dm_control instance
1185  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
1186  CHKERR DMSetType(dm_sub_force, dm_name);
1187  // set dm_sub_force data structure which created mofem data structures
1188  CHKERR DMMoFEMCreateSubDM(dm_sub_force, dm_control, "FORCES_PROB");
1189  CHKERR DMMoFEMSetSquareProblem(dm_sub_force, PETSC_FALSE);
1190  CHKERR DMMoFEMAddSubFieldRow(dm_sub_force, "RHO");
1191  CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "U");
1192  CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "UPSILON");
1193  // add elements to dm_sub_force
1194  CHKERR DMMoFEMAddElement(dm_sub_force, "B");
1195  CHKERR DMSetUp(dm_sub_force);
1196 
1197  Mat UB, UPSILONB;
1198  CHKERR DMCreateMatrix(dm_sub_force, &UB);
1199  CHKERR MatZeroEntries(UB);
1200  // note be will be transposed in place latter on
1201  CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
1202  CHKERR MatZeroEntries(UPSILONB);
1203  {
1204  CellEngineering::FaceElement face_b_matrices(m_field);
1205  face_b_matrices.getOpPtrVector().push_back(
1206  new OpCalculateInvJacForFace(inv_jac));
1207  if (is_curl) {
1208  face_b_matrices.getOpPtrVector().push_back(
1209  new OpSetInvJacH1ForFace(inv_jac));
1210  face_b_matrices.getOpPtrVector().push_back(
1211  new OpSetInvJacHcurlFace(inv_jac));
1212  face_b_matrices.getOpPtrVector().push_back(new OpCellCurlB(UB, "U"));
1213  face_b_matrices.getOpPtrVector().push_back(
1214  new OpCellCurlB(UPSILONB, "UPSILON"));
1215  } else {
1216  face_b_matrices.getOpPtrVector().push_back(
1217  new OpSetInvJacH1ForFace(inv_jac));
1218  face_b_matrices.getOpPtrVector().push_back(
1219  new OpCellPotentialB(UB, "U"));
1220  face_b_matrices.getOpPtrVector().push_back(
1221  new OpCellPotentialB(UPSILONB, "UPSILON"));
1222  }
1223  CHKERR DMoFEMLoopFiniteElements(dm_sub_force, "B", &face_b_matrices);
1224  }
1225  CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
1226  CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
1227  CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
1228  CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
1229 
1230  const Problem *problem_ptr;
1231  CHKERR m_field.get_problem("FORCES_PROB", &problem_ptr);
1232 
1233  // Zero rows, force field is given by gradients of potential field, so one
1234  // of the values has to be fixed like for rigid body motion.
1235  if (is_curl == PETSC_FALSE) {
1236  int nb_dofs_to_fix = 0;
1237  int index_to_fix = 0;
1238  if (!vertex_to_fix.empty()) {
1239  boost::shared_ptr<NumeredDofEntity> dof_ptr;
1240  CHKERR problem_ptr->getDofByNameEntAndEntDofIdx(
1241  m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
1242  dof_ptr);
1243  if (dof_ptr) {
1244  if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1245  nb_dofs_to_fix = 1;
1246  index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1247  cerr << *dof_ptr << endl;
1248  }
1249  }
1250  }
1251  CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
1252  PETSC_NULL);
1253  CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
1254  PETSC_NULL, PETSC_NULL);
1255  } else {
1256  std::vector<int> dofs_to_fix;
1257  for (auto p_eit = edges_to_fix.pair_begin();
1258  p_eit != edges_to_fix.pair_end(); ++p_eit) {
1259  auto bit_number = m_field.get_field_bit_number("RHO");
1260  auto row_dofs = problem_ptr->numeredRowDofsPtr;
1261  auto lo = row_dofs->lower_bound(
1262  FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
1263  auto hi =
1264  row_dofs->lower_bound(FieldEntity::getHiLocalEntityBitNumber(
1265  bit_number, p_eit->second));
1266  for (; lo != hi; ++lo)
1267  dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1268  }
1269  CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
1270  PETSC_NULL, PETSC_NULL);
1271  CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1272  0, PETSC_NULL, PETSC_NULL);
1273  }
1274 
1275  Mat UBT;
1276  CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
1277  CHKERR MatDestroy(&UB);
1278 
1279  // Matrix View
1280  if (debug) {
1281  cerr << "UBT" << endl;
1282  MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
1283  std::string wait;
1284  std::cin >> wait;
1285  }
1286 
1287  boost::shared_ptr<Problem::SubProblemData> sub_data =
1288  problem_ptr->getSubData();
1289  // CHKERR sub_data->getColIs(&nested_is_rows[0]);
1290  // CHKERR sub_data->getRowIs(&nested_is_cols[1]);
1291  nested_matrices(0, 1) = UBT;
1292 
1293  if (debug) {
1294  cerr << "UPSILONB" << endl;
1295  MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
1296  std::string wait;
1297  std::cin >> wait;
1298  }
1299 
1300  // CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1301  // CHKERR sub_data->getColIs(&nested_is_cols[0]);
1302  nested_matrices(1, 0) = UPSILONB;
1303 
1304  CHKERR DMDestroy(&dm_sub_force);
1305  }
1306 
1307  Mat SubA;
1308  CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
1309  &sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
1310  &SubA);
1311  nested_matrices(0, 0) = SubA;
1312 
1313  CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
1314  CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
1315 
1316  if (debug) {
1317  cerr << "Nested SubA" << endl;
1318  MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
1319  }
1320 
1321  Mat A;
1322  CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
1323  &nested_is_cols[0], &nested_matrices(0, 0), &A);
1324 
1325  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1326  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1327 
1328  if (debug) {
1329  cerr << "Nested A" << endl;
1330  MatView(A, PETSC_VIEWER_STDOUT_WORLD);
1331  }
1332 
1333  Vec D, F;
1334  CHKERR DMCreateGlobalVector(dm_control, &D);
1335  CHKERR DMCreateGlobalVector(dm_control, &F);
1336 
1337  // Asseble the right hand side vector
1338  {
1339  CellEngineering::FaceElement face_element(m_field);
1340  face_element.getOpPtrVector().push_back(new OpGetDispX(common_data));
1341  face_element.getOpPtrVector().push_back(new OpGetDispY(common_data));
1342  face_element.getOpPtrVector().push_back(
1343  new OpCell_g(F, eps_u, common_data));
1344  CHKERR DMoFEMLoopFiniteElements(dm_control, "DISPLACEMENTS_PENALTY",
1345  &face_element);
1346  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
1347  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
1348  CHKERR VecAssemblyBegin(F);
1349  CHKERR VecAssemblyEnd(F);
1350  }
1351 
1352  KSP solver;
1353  // KSP solver;
1354  {
1355  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
1356  CHKERR KSPSetDM(solver, dm_control);
1357  CHKERR KSPSetFromOptions(solver);
1358  CHKERR KSPSetOperators(solver, A, A);
1359  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
1360  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
1361  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
1362  PC pc;
1363  CHKERR KSPGetPC(solver, &pc);
1364  CHKERR PCSetType(pc, PCFIELDSPLIT);
1365  PetscBool is_pcfs = PETSC_FALSE;
1366  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1367  if (is_pcfs) {
1368  CHKERR PCSetOperators(pc, A, A);
1369  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
1370  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
1371  CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1372  CHKERR PCSetUp(pc);
1373  KSP *sub_ksp;
1374  PetscInt n;
1375  CHKERR PCFieldSplitGetSubKSP(pc, &n, &sub_ksp);
1376  {
1377  PC sub_pc_0;
1378  CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
1379  CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
1380  CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
1381  CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
1382  CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
1383  CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
1384  // CHKERR
1385  // PCFieldSplitSetSchurFactType(sub_pc_0,PC_FIELDSPLIT_SCHUR_FACT_LOWER);
1386  // CHKERR PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_SCHUR);
1387  CHKERR PCSetUp(sub_pc_0);
1388  }
1389  } else {
1390  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1391  "Only works with pre-conditioner PCFIELDSPLIT");
1392  }
1393  CHKERR KSPSetUp(solver);
1394  }
1395 
1396  // Solve system of equations
1397  CHKERR KSPSolve(solver, F, D);
1398 
1399  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1400  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1401  CHKERR DMoFEMMeshToLocalVector(dm_control, D, INSERT_VALUES,
1402  SCATTER_REVERSE);
1403 
1404  if (debug) {
1405  CHKERR VecView(D, PETSC_VIEWER_DRAW_WORLD);
1406  std::string wait;
1407  std::cin >> wait;
1408  }
1409 
1410  // Clean sub matrices and sub indices
1411  for (int i = 0; i != 2; i++) {
1412  if (sub_nested_is_rows[i]) {
1413  CHKERR ISDestroy(&sub_nested_is_rows[i]);
1414  }
1415  if (sub_nested_is_cols[i]) {
1416  CHKERR ISDestroy(&sub_nested_is_cols[i]);
1417  }
1418  for (int j = 0; j != 2; j++) {
1419  if (sub_nested_matrices(i, j)) {
1420  CHKERR MatDestroy(&sub_nested_matrices(i, j));
1421  }
1422  }
1423  }
1424  for (int i = 0; i != 2; i++) {
1425  if (nested_is_rows[i]) {
1426  CHKERR ISDestroy(&nested_is_rows[i]);
1427  }
1428  if (nested_is_cols[i]) {
1429  CHKERR ISDestroy(&nested_is_cols[i]);
1430  }
1431  for (int j = 0; j != 2; j++) {
1432  if (nested_matrices(i, j)) {
1433  CHKERR MatDestroy(&nested_matrices(i, j));
1434  }
1435  }
1436  }
1437 
1438  CHKERR MatDestroy(&SubA);
1439  CHKERR MatDestroy(&A);
1440  CHKERR VecDestroy(&D);
1441  CHKERR VecDestroy(&F);
1442 
1443  CHKERR DMDestroy(&dm_sub_volume_control);
1444 
1445  PostProcVolumeOnRefinedMesh post_proc(m_field);
1446  {
1448  CHKERR post_proc.addFieldValuesPostProc("U");
1449  CHKERR post_proc.addFieldValuesGradientPostProc("U");
1450  // add postpocessing for sresses
1451  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
1452  m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "U",
1453  post_proc.commonData, &elastic.setOfBlocks));
1454  CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC", &post_proc);
1455  CHKERR post_proc.writeFile("out.h5m");
1456  elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
1457  elastic.getLoopFeEnergy().eNergy = 0;
1458  CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC",
1459  &elastic.getLoopFeEnergy());
1460  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
1461  elastic.getLoopFeEnergy().eNergy);
1462  }
1463 
1464  {
1465  PostProcFaceOnRefinedMesh post_proc_face(m_field);
1466  CHKERR post_proc_face.generateReferenceElementMesh();
1467  post_proc_face.getOpPtrVector().push_back(
1468  new OpCalculateInvJacForFace(inv_jac));
1469  if (is_curl) {
1470  post_proc_face.getOpPtrVector().push_back(
1471  new OpSetInvJacHcurlFace(inv_jac));
1472  post_proc_face.getOpPtrVector().push_back(
1473  new OpVirtualCurlRho("RHO", common_data));
1474  } else {
1475  post_proc_face.getOpPtrVector().push_back(
1476  new OpSetInvJacH1ForFace(inv_jac));
1477  post_proc_face.getOpPtrVector().push_back(
1478  new OpVirtualPotentialRho("RHO", common_data));
1479  }
1480  post_proc_face.getOpPtrVector().push_back(
1481  new PostProcTraction(m_field, post_proc_face.postProcMesh,
1482  post_proc_face.mapGaussPts, common_data));
1483  CHKERR DMoFEMLoopFiniteElements(dm_control, "D", &post_proc_face);
1484  CHKERR post_proc_face.writeFile("out_tractions.h5m");
1485  }
1486 
1487  CHKERR DMDestroy(&dm_control);
1488 
1489  }
1490  CATCH_ERRORS;
1491 
1493  return 0;
1494 }
CellEngineering::BlockOptionData::yOung
double yOung
Definition: cell_forces.cpp:416
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
SIDESET
@ SIDESET
Definition: definitions.h:160
debug
static int debug
Definition: cell_forces.cpp:430
CellEngineering::PostProcTraction
Shave results on mesh tags for post-processing.
Definition: CellForces.hpp:600
CellEngineering
Definition: cell_forces.cpp:412
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::OpSetInvJacHcurlFace
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
Definition: UserDataOperators.hpp:2978
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress
Operator performs automatic differentiation.
Definition: NonLinearElasticElement.hpp:370
MoFEM::UnknownInterface::query_interface
virtual MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
CellEngineering::FaceElement
Definition: CellForces.hpp:48
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::FieldEntity::getLoLocalEntityBitNumber
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition: FieldEntsMultiIndices.hpp:247
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition: LoopMethods.hpp:124
CellEngineering::OpCellCurlB
Calculate and assemble Z matrix.
Definition: CellForces.hpp:392
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DMoFEMMeshToLocalVector
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:523
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3145
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
help
static char help[]
Definition: cell_forces.cpp:409
CellEngineering::OpVirtualCurlRho
Post-process tractions.
Definition: CellForces.hpp:528
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:304
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3172
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::MeshsetsManager::getEntitiesByDimension
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimension
Definition: MeshsetsManager.cpp:669
main
int main(int argc, char *argv[])
Definition: cell_forces.cpp:432
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
Hooke.hpp
Implementation of linear elastic material.
MoFEM::DMMoFEMCreateSubDM
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:215
PostProcTemplateVolumeOnRefinedMesh::commonData
CommonData commonData
Definition: PostProcOnRefMesh.hpp:287
CellEngineering::OpCellCurlD
Calculate and assemble Z matrix.
Definition: CellForces.hpp:315
CellEngineering::OpCellPotentialD
Calculate and assemble D matrix.
Definition: CellForces.hpp:73
MoFEM::Problem::getSubData
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
Definition: ProblemsMultiIndices.hpp:119
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:556
PostProcTemplateOnRefineMesh::writeFile
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"
Definition: PostProcOnRefMesh.hpp:253
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::TSMethod::CTX_TSNONE
@ CTX_TSNONE
Definition: LoopMethods.hpp:145
CellEngineering::BlockOptionData
Definition: cell_forces.cpp:414
MoFEM::DeprecatedCoreInterface::synchronise_entities
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:41
CellEngineering::FatPrism
Definition: CellForces.hpp:21
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Definition: PostProcOnRefMesh.hpp:195
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
NonlinearElasticElement::addElement
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)
Definition: NonLinearElasticElement.cpp:1120
CellEngineering::OpCellS
Calculate and assemble S matrix.
Definition: CellForces.hpp:201
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
CellEngineering::OpGetDispX
Definition: CellForces.hpp:57
MoFEM::DMMoFEMCreateMoFEM
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:114
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:325
NonlinearElasticElement::setOperators
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.
Definition: NonLinearElasticElement.cpp:1153
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
convert.n
n
Definition: convert.py:82
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition: ProblemsMultiIndices.hpp:73
NonlinearElasticElement::commonData
CommonData commonData
Definition: NonLinearElasticElement.hpp:121
MoFEM::OpSetInvJacH1ForFace
Definition: UserDataOperators.hpp:2918
Range
CellEngineering::BlockOptionData::BlockOptionData
BlockOptionData()
Definition: cell_forces.cpp:418
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
MoFEM::CoreTmp< 0 >::Initialize
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
CellEngineering::BlockOptionData::oRder
int oRder
Definition: cell_forces.cpp:415
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
MoFEM::DeprecatedCoreInterface::seed_ref_level
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
Definition: DeprecatedCoreInterface.cpp:24
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:290
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
std
Definition: enable_if.hpp:5
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
CellEngineering::OpCellPotentialB
Calculate and assemble B matrix.
Definition: CellForces.hpp:140
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
CellEngineering::CommonData
Definition: CellForces.hpp:31
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
CellEngineering::OpGetDispY
Definition: CellForces.hpp:65
MoFEM::DeprecatedCoreInterface::synchronise_field_entities
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:47
MoFEM::FieldEntity::getHiLocalEntityBitNumber
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition: FieldEntsMultiIndices.hpp:258
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
MoFEM::CoreInterface::set_field_order
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.
MoFEM::Problem::getDofByNameEntAndEntDofIdx
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
Definition: ProblemsMultiIndices.cpp:132
CellEngineering::BlockOptionData::pOisson
double pOisson
Definition: cell_forces.cpp:417
MoFEM::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:360
CellEngineering::OpCell_g
Calculate and assemble g vector.
Definition: CellForces.hpp:267
CellEngineering::OpVirtualPotentialRho
Post-process tractions.
Definition: CellForces.hpp:460
MoFEM::DMoFEMLoopFiniteElements
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:586
MoFEM::CoreInterface::add_field
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.
PostProcFaceOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: PostProcOnRefMesh.cpp:539
NonlinearElasticElement::setBlocks
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double >> materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble >> materialAdoublePtr)
Definition: NonLinearElasticElement.cpp:1086
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
F
@ F
Definition: free_surface.cpp:394
Hooke
Hook equation.
Definition: Hooke.hpp:19
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
CellForces.hpp
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition: PostProcOnRefMesh.hpp:153