v0.14.0
stability.cpp
Go to the documentation of this file.
1 /** \file stability.cpp
2  * \ingroup nonlinear_elastic_elem
3  *
4  * Solves stability problem. Currently uses 3d tetrahedral elements.
5  */
6 
7 
8 
9 #include <BasicFiniteElements.hpp>
10 using namespace MoFEM;
11 #include <Hooke.hpp>
12 
13 #undef EPS
14 #include <slepceps.h>
15 
16 #include <SurfacePressureComplexForLazy.hpp>
17 
18 static char help[] = "...\n\n";
19 
20 template <typename TYPE>
21 struct MyMat_double
22  : public NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<
23  TYPE> {
24 
25  bool doAotherwiseB;
26  MyMat_double() : doAotherwiseB(true){};
27 
28  MatrixDouble D_lambda, D_mu, D;
29  VectorBoundedArray<TYPE, 6> sTrain, sTrain0, sTress;
30 
31  FTensor::Index<'i', 3> i;
32  FTensor::Index<'j', 3> j;
33  FTensor::Index<'k', 3> k;
34 
35  virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(
36  const NonlinearElasticElement::BlockData block_data,
37  boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
38  MoFEMFunctionBegin;
39  double lambda = LAMBDA(block_data.E, block_data.PoissonRatio);
40  double mu = MU(block_data.E, block_data.PoissonRatio);
41  if (D_lambda.size1() == 0) {
42  D_lambda.resize(6, 6);
43  D_lambda.clear();
44  for (int rr = 0; rr < 3; rr++) {
45  for (int cc = 0; cc < 3; cc++) {
46  D_lambda(rr, cc) = 1;
47  }
48  }
49  }
50  if (D_mu.size1() == 0) {
51  D_mu.resize(6, 6);
52  D_mu.clear();
53  for (int rr = 0; rr < 6; rr++) {
54  D_mu(rr, rr) = rr < 3 ? 2 : 1;
55  }
56  }
57  D.resize(6, 6);
58  noalias(D) = lambda * D_lambda + mu * D_mu;
59 
60  if (doAotherwiseB) {
61  sTrain.resize(6);
62  sTress.resize(6);
63 
64  sTrain[0] = this->F(0, 0) - 1;
65  sTrain[1] = this->F(1, 1) - 1;
66  sTrain[2] = this->F(2, 2) - 1;
67  sTrain[3] = this->F(0, 1) + this->F(1, 0);
68  sTrain[4] = this->F(1, 2) + this->F(2, 1);
69  sTrain[5] = this->F(0, 2) + this->F(2, 0);
70 
71  sTress.clear();
72  for (int ii = 0; ii != 6; ++ii)
73  for (int jj = 0; jj != 6; ++jj)
74  sTress(ii) += D(ii, jj) * sTrain(jj);
75 
76  this->P(0, 0) = sTress[0];
77  this->P(1, 1) = sTress[1];
78  this->P(2, 2) = sTress[2];
79  this->P(0, 1) = this->P(1, 0) = sTress[3];
80  this->P(1, 2) = this->P(2, 1) = sTress[4];
81  this->P(0, 2) = this->P(2, 0) = sTress[5];
82 
83  } else {
84 
85  sTrain0.resize(6, 0);
86 
87  sTress.clear();
88  for (int ii = 0; ii != 6; ++ii)
89  for (int jj = 0; jj != 6; ++jj)
90  sTress(ii) += D(ii, jj) * sTrain(jj);
91 
92  // noalias(sTress) = prod(D, sTrain0);
93  this->sigmaCauchy(0, 0) = sTress[0];
94  this->sigmaCauchy(1, 1) = sTress[1];
95  this->sigmaCauchy(2, 2) = sTress[2];
96  this->sigmaCauchy(0, 1) = this->sigmaCauchy(1, 0) = sTress[3];
97  this->sigmaCauchy(1, 2) = this->sigmaCauchy(2, 1) = sTress[4];
98  this->sigmaCauchy(0, 2) = this->sigmaCauchy(2, 0) = sTress[5];
99 
100  auto J = determinantTensor3by3(this->F);
101  CHKERR invertTensor3by3(this->F, J, this->invF);
102 
103  this->t_P(i, j) = this->t_sigmaCauchy(i, k) * this->t_invF(j, k);
104  this->t_P(i, j) *= J;
105  }
106 
107  MoFEMFunctionReturn(0);
108  }
109 };
110 
111 template <typename TYPE> struct MyMat : public MyMat_double<TYPE> {
112 
113  int nbActiveVariables0;
114 
115  virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables) {
116  MoFEMFunctionBeginHot;
117 
118  try {
119 
120  this->sTrain0.resize(6);
121  MatrixDouble &G0 = (this->commonDataPtr->gradAtGaussPts["D0"][this->gG]);
122  this->sTrain0[0] <<= G0(0, 0);
123  this->sTrain0[1] <<= G0(1, 1);
124  this->sTrain0[2] <<= G0(2, 2);
125  this->sTrain0[3] <<= (G0(1, 0) + G0(0, 1));
126  this->sTrain0[4] <<= (G0(2, 1) + G0(1, 2));
127  this->sTrain0[5] <<= (G0(2, 0) + G0(0, 2));
128  nbActiveVariables0 = nb_active_variables;
129  nb_active_variables += 6;
130 
131  } catch (const std::exception &ex) {
132  std::ostringstream ss;
133  ss << "throw in method: " << ex.what() << std::endl;
134  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
135  }
136 
137  MoFEMFunctionReturnHot(0);
138  }
139 
140  virtual MoFEMErrorCode
141  setUserActiveVariables(VectorDouble &active_variables) {
142  MoFEMFunctionBeginHot;
143 
144  try {
145 
146  int shift = nbActiveVariables0; // is a number of elements in F
147  MatrixDouble &G0 = (this->commonDataPtr->gradAtGaussPts["D0"][this->gG]);
148  active_variables[shift + 0] = G0(0, 0);
149  active_variables[shift + 1] = G0(1, 1);
150  active_variables[shift + 2] = G0(2, 2);
151  active_variables[shift + 3] = G0(0, 1) + G0(1, 0);
152  active_variables[shift + 4] = G0(1, 2) + G0(2, 1);
153  active_variables[shift + 5] = G0(0, 2) + G0(2, 0);
154 
155  } catch (const std::exception &ex) {
156  std::ostringstream ss;
157  ss << "throw in method: " << ex.what() << std::endl;
158  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
159  }
160 
161  MoFEMFunctionReturnHot(0);
162  }
163 };
164 
165 int main(int argc, char *argv[]) {
166 
167  // PetscInitialize(&argc,&argv,(char *)0,help);
168  const string default_options = "-ksp_type fgmres \n"
169  "-pc_type lu \n"
170  "-pc_factor_mat_solver_type mumps \n"
171  "-mat_mumps_icntl_20 0 \n"
172  "-ksp_atol 1e-10 \n"
173  "-ksp_rtol 1e-10 \n"
174  "-snes_monitor \n"
175  "-snes_type newtonls \n"
176  "-snes_linesearch_type basic \n"
177  "-snes_max_it 100 \n"
178  "-snes_atol 1e-7 \n"
179  "-snes_rtol 1e-7 \n"
180  "-ts_monitor \n"
181  "-ts_type alpha \n";
182 
183  string param_file = "param_file.petsc";
184  if (!static_cast<bool>(ifstream(param_file))) {
185  std::ofstream file(param_file.c_str(), std::ios::ate);
186  if (file.is_open()) {
187  file << default_options;
188  file.close();
189  }
190  }
191 
192  SlepcInitialize(&argc, &argv, param_file.c_str(), help);
193  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
194 
195  try {
196 
197  moab::Core mb_instance;
198  moab::Interface &moab = mb_instance;
199  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
200  auto moab_comm_wrap =
201  boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
202  if (pcomm == NULL)
203  pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
204 
205  PetscBool flg = PETSC_TRUE;
206  char mesh_file_name[255];
207  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
208  mesh_file_name, 255, &flg);
209  if (flg != PETSC_TRUE) {
210  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
211  }
212 
213  // use this if your mesh is partotioned and you run code on parts,
214  // you can solve very big problems
215  PetscBool is_partitioned = PETSC_FALSE;
216  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_partitioned",
217  &is_partitioned, &flg);
218 
219  if (is_partitioned == PETSC_TRUE) {
220  // Read mesh to MOAB
221  const char *option;
222  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
223  "PARALLEL_PARTITION;";
224  CHKERR moab.load_file(mesh_file_name, 0, option);
225  CHKERR pcomm->resolve_shared_ents(0, 3, 0);
226  CHKERR pcomm->resolve_shared_ents(0, 3, 1);
227  CHKERR pcomm->resolve_shared_ents(0, 3, 2);
228  } else {
229  const char *option;
230  option = "";
231  CHKERR moab.load_file(mesh_file_name, 0, option);
232  }
233 
234  MoFEM::Core core(moab);
235  MoFEM::Interface &m_field = core;
236 
237  Range CubitSIDESETs_meshsets;
238  CHKERR m_field.getInterface<MeshsetsManager>()->getMeshsetsByType(
239  SIDESET, CubitSIDESETs_meshsets);
240 
241  // ref meshset ref level 0
242  BitRefLevel bit_level0;
243  bit_level0.set(0);
244  EntityHandle meshset_level0;
245  CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
246  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
247  0, 3, bit_level0);
248  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
249  bit_level0, BitRefLevel().set(), meshset_level0);
250 
251  // Fields
252  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
253  3, MB_TAG_SPARSE, MF_ZERO);
254  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
255  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
256  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
257  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
258  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
259 
260  bool check_if_spatial_field_exist = m_field.check_field("SPATIAL_POSITION");
261  CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
262  MB_TAG_SPARSE, MF_ZERO);
263  CHKERR m_field.add_field("EIGEN_VECTOR", H1, AINSWORTH_LEGENDRE_BASE, 3,
264  MB_TAG_SPARSE, MF_ZERO);
265  CHKERR m_field.add_field("D0", H1, AINSWORTH_LEGENDRE_BASE, 3,
266  MB_TAG_SPARSE, MF_ZERO);
267 
268  // add entitities (by tets) to the field
269  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
270  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_VECTOR");
271  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "D0");
272 
273  boost::shared_ptr<Hooke<double>> mat_double =
274  boost::make_shared<Hooke<double>>();
275  boost::shared_ptr<MyMat<adouble>> mat_adouble =
276  boost::make_shared<MyMat<adouble>>();
277 
278  NonlinearElasticElement elastic(m_field, 2);
279  CHKERR elastic.setBlocks(mat_double, mat_adouble);
280  CHKERR elastic.addElement("ELASTIC", "SPATIAL_POSITION");
281  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC",
282  "EIGEN_VECTOR");
283  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "D0");
284 
285  elastic.feRhs.getOpPtrVector().push_back(
286  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
287  "D0", elastic.commonData));
288  elastic.feLhs.getOpPtrVector().push_back(
289  new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
290  "D0", elastic.commonData));
291  CHKERR elastic.setOperators("SPATIAL_POSITION");
292 
293  // define problems
294  CHKERR m_field.add_problem("ELASTIC_MECHANICS", MF_ZERO);
295  // set finite elements for problems
296  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
297  "ELASTIC");
298  // set refinement level for problem
299  CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
300  bit_level0);
301 
302  // set app. order
303 
304  PetscInt disp_order;
305  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &disp_order,
306  &flg);
307  if (flg != PETSC_TRUE) {
308  disp_order = 1;
309  }
310 
311  CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", disp_order);
312  CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", disp_order);
313  CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", disp_order);
314  CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
315  CHKERR m_field.set_field_order(0, MBTET, "EIGEN_VECTOR", disp_order);
316  CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_VECTOR", disp_order);
317  CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_VECTOR", disp_order);
318  CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_VECTOR", 1);
319  CHKERR m_field.set_field_order(0, MBTET, "D0", disp_order);
320  CHKERR m_field.set_field_order(0, MBTRI, "D0", disp_order);
321  CHKERR m_field.set_field_order(0, MBEDGE, "D0", disp_order);
322  CHKERR m_field.set_field_order(0, MBVERTEX, "D0", 1);
323 
324  CHKERR m_field.add_finite_element("NEUAMNN_FE", MF_ZERO);
325  CHKERR m_field.modify_finite_element_add_field_row("NEUAMNN_FE",
326  "SPATIAL_POSITION");
327  CHKERR m_field.modify_finite_element_add_field_col("NEUAMNN_FE",
328  "SPATIAL_POSITION");
329  CHKERR m_field.modify_finite_element_add_field_data("NEUAMNN_FE",
330  "SPATIAL_POSITION");
331  CHKERR m_field.modify_finite_element_add_field_data("NEUAMNN_FE",
332  "MESH_NODE_POSITIONS");
333  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
334  "NEUAMNN_FE");
335  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
336  it)) {
337  Range tris;
338  CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, true);
339  CHKERR m_field.add_ents_to_finite_element_by_type(tris, MBTRI,
340  "NEUAMNN_FE");
341  }
342  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
343  m_field, SIDESET | PRESSURESET, it)) {
344  Range tris;
345  CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, true);
346  CHKERR m_field.add_ents_to_finite_element_by_type(tris, MBTRI,
347  "NEUAMNN_FE");
348  }
349  // add nodal force element
350  CHKERR MetaNodalForces::addElement(m_field, "SPATIAL_POSITION");
351  CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
352  "FORCE_FE");
353 
354  // build field
355  CHKERR m_field.build_fields();
356  // 10 node tets
357  if (!check_if_spatial_field_exist) {
358  Projection10NodeCoordsOnField ent_method_material(m_field,
359  "MESH_NODE_POSITIONS");
360  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
361  Projection10NodeCoordsOnField ent_method_spatial(m_field,
362  "SPATIAL_POSITION");
363  CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method_spatial);
364  // CHKERR m_field.set_field(0,MBTRI,"SPATIAL_POSITION");
365  // CHKERR m_field.set_field(0,MBTET,"SPATIAL_POSITION");
366  // CHKERR m_field.field_axpy(1,"SPATIAL_POSITION","D0",true);
367  }
368 
369  // build finite elemnts
370  CHKERR m_field.build_finite_elements();
371  // build adjacencies
372  CHKERR m_field.build_adjacencies(bit_level0);
373 
374  // build database
375  ProblemsManager *prb_mng_ptr;
376  CHKERR m_field.getInterface(prb_mng_ptr);
377  if (is_partitioned) {
378  CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("ELASTIC_MECHANICS",
379  true);
380  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS", true, 0,
381  pcomm->size(), 1);
382  } else {
383  CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
384  CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
385  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS");
386  }
387  CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
388 
389  // create matrices
390  Vec F;
391  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
392  "ELASTIC_MECHANICS", ROW, &F);
393  Vec D;
394  CHKERR VecDuplicate(F, &D);
395  Mat Aij;
396  CHKERR m_field.getInterface<MatrixManager>()
397  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
398  &Aij);
399 
400  // surface forces
401  NeumannForcesSurfaceComplexForLazy neumann_forces(m_field, Aij, F);
402  NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE &neumann =
403  neumann_forces.getLoopSpatialFe();
404  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, NODESET | FORCESET,
405  it)) {
406  CHKERR neumann.addForce(it->getMeshsetId());
407  }
408  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
409  m_field, SIDESET | PRESSURESET, it)) {
410  CHKERR neumann.addPressure(it->getMeshsetId());
411  }
412  DirichletSpatialPositionsBc my_Dirichlet_bc(m_field, "SPATIAL_POSITION",
413  Aij, D, F);
414 
415  CHKERR VecZeroEntries(F);
416  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
417  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
418  CHKERR MatZeroEntries(Aij);
419 
420  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
421  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
422  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
423  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
424 
425  // F Vector
426  // preproc
427  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
428  my_Dirichlet_bc.snes_x = D;
429  my_Dirichlet_bc.snes_f = F;
430  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
431  my_Dirichlet_bc);
432  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
433  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
434  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
435  "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
436  // elem loops
437  // noadl forces
438  boost::ptr_map<std::string, NodalForce> nodal_forces;
439  string fe_name_str = "FORCE_FE";
440  nodal_forces.insert(fe_name_str, new NodalForce(m_field));
441  CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
442  "SPATIAL_POSITION");
443  boost::ptr_map<std::string, NodalForce>::iterator fit =
444  nodal_forces.begin();
445  for (; fit != nodal_forces.end(); fit++) {
446  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", fit->first,
447  fit->second->getLoopFe());
448  }
449  // surface forces
450  neumann.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
451  neumann.snes_x = D;
452  neumann.snes_f = F;
453  m_field.loop_finite_elements("ELASTIC_MECHANICS", "NEUAMNN_FE", neumann);
454  // stiffnes
455  elastic.getLoopFeRhs().snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
456  elastic.getLoopFeRhs().snes_x = D;
457  elastic.getLoopFeRhs().snes_f = F;
458  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
459  elastic.getLoopFeRhs());
460  // postproc
461  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
462  my_Dirichlet_bc);
463 
464  // Aij Matrix
465  // preproc
466  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
467  my_Dirichlet_bc.snes_B = Aij;
468  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
469  my_Dirichlet_bc);
470  // surface forces
471  // neumann.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
472  // neumann.snes_B = Aij;
473  // CHKERR
474  // m_field.loop_finite_elements("ELASTIC_MECHANICS","NEUAMNN_FE",neumann);
475  // stiffnes
476  elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
477  elastic.getLoopFeLhs().snes_B = Aij;
478  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
479  elastic.getLoopFeLhs());
480  // postproc
481  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
482  my_Dirichlet_bc);
483 
484  CHKERR VecAssemblyBegin(F);
485  CHKERR VecAssemblyEnd(F);
486  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
487  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
488 
489  CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
490  CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
491 
492  // Matrix View
493  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
494  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
495  // std::string wait;
496  // std::cin >> wait;
497 
498  // Solver
499  KSP solver;
500  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
501  CHKERR KSPSetOperators(solver, Aij, Aij);
502  CHKERR KSPSetFromOptions(solver);
503 
504  CHKERR KSPSetUp(solver);
505 
506  CHKERR VecZeroEntries(D);
507  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
508  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
509 
510  CHKERR KSPSolve(solver, F, D);
511  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
512  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
513 
514  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
515  "ELASTIC_MECHANICS", "SPATIAL_POSITION", "D0", COL, D, INSERT_VALUES,
516  SCATTER_REVERSE);
517 
518  Mat Bij;
519  CHKERR MatDuplicate(Aij, MAT_SHARE_NONZERO_PATTERN, &Bij);
520  // CHKERR MatZeroEntries(Aij);
521  CHKERR MatZeroEntries(Bij);
522 
523  /*//Aij Matrix
524  //preproc
525  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
526  my_Dirichlet_bc.snes_B = Aij;
527  CHKERR
528  m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",my_Dirichlet_bc);
529  //stiffnes
530  elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
531  elastic.getLoopFeLhs().snes_B = Aij;
532  CHKERR
533  m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs());
534  //postproc
535  CHKERR
536  m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",my_Dirichlet_bc);
537  */
538 
539  // Bij Matrix
540  mat_adouble->doAotherwiseB = false;
541  // preproc
542  my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
543  my_Dirichlet_bc.snes_B = Bij;
544  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
545  my_Dirichlet_bc);
546  // surface forces
547  neumann.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
548  neumann.snes_B = Bij;
549  PetscBool is_conservative = PETSC_FALSE;
550  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_conservative",
551  &is_conservative, &flg);
552  if (is_conservative) {
553  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "NEUAMNN_FE",
554  neumann);
555  }
556  // stiffnes
557  elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
558  elastic.getLoopFeLhs().snes_B = Bij;
559  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
560  elastic.getLoopFeLhs());
561  // postproc
562  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
563  my_Dirichlet_bc);
564 
565  CHKERR MatSetOption(Bij, MAT_SPD, PETSC_TRUE);
566  CHKERR MatAssemblyBegin(Bij, MAT_FINAL_ASSEMBLY);
567  CHKERR MatAssemblyEnd(Bij, MAT_FINAL_ASSEMBLY);
568 
569  // Matrix View
570  // MatView(Bij,PETSC_VIEWER_STDOUT_WORLD);
571  // MatView(Bij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
572  // std::string wait;
573  // std::cin >> wait;
574 
575  EPS eps;
576  ST st;
577  EPSType type;
578  PetscReal tol;
579  PetscInt nev, maxit, its;
580 
581  /*
582  Create eigensolver context
583  */
584  CHKERR EPSCreate(PETSC_COMM_WORLD, &eps);
585  /*
586  Set operators. In this case, it is a generalized eigenvalue problem
587  */
588  CHKERR EPSSetOperators(eps, Bij, Aij);
589  /*
590  Set solver parameters at runtime
591  */
592  CHKERR EPSSetFromOptions(eps);
593  /*
594  Optional: Get some information from the solver and display it
595  */
596  CHKERR EPSSolve(eps);
597 
598  /*
599  Optional: Get some information from the solver and display it
600  */
601  CHKERR EPSGetIterationNumber(eps, &its);
602  PetscPrintf(PETSC_COMM_WORLD, " Number of iterations of the method: %D\n",
603  its);
604  CHKERR EPSGetST(eps, &st);
605  // CHKERR STGetOperationCounters(st,NULL,&lits);
606  // PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method:
607  // %D\n",lits);
608  CHKERR EPSGetType(eps, &type);
609  PetscPrintf(PETSC_COMM_WORLD, " Solution method: %s\n", type);
610  CHKERR EPSGetDimensions(eps, &nev, NULL, NULL);
611  PetscPrintf(PETSC_COMM_WORLD, " Number of requested eigenvalues: %D\n",
612  nev);
613  CHKERR EPSGetTolerances(eps, &tol, &maxit);
614  PetscPrintf(PETSC_COMM_WORLD, " Stopping condition: tol=%.4g, maxit=%D\n",
615  (double)tol, maxit);
616 
617  // get solutions
618  PostProcVolumeOnRefinedMesh post_proc(m_field);
619  CHKERR post_proc.generateReferenceElementMesh();
620  CHKERR post_proc.addFieldValuesGradientPostProc("SPATIAL_POSITION");
621  CHKERR post_proc.addFieldValuesPostProc("SPATIAL_POSITION");
622  CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
623  CHKERR post_proc.addFieldValuesPostProc("EIGEN_VECTOR");
624  CHKERR post_proc.addFieldValuesGradientPostProc("EIGEN_VECTOR");
625  CHKERR post_proc.addFieldValuesPostProc("D0");
626  CHKERR post_proc.addFieldValuesGradientPostProc("D0");
627  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
628  post_proc);
629  CHKERR post_proc.writeFile("out.h5m");
630 
631  PetscScalar eigr, eigi, nrm2r;
632  for (int nn = 0; nn < nev; nn++) {
633  CHKERR EPSGetEigenpair(eps, nn, &eigr, &eigi, D, PETSC_NULL);
634  CHKERR VecNorm(D, NORM_2, &nrm2r);
635  PetscPrintf(
636  PETSC_COMM_WORLD,
637  " ncov = %D eigr = %.4g eigi = %.4g (inv eigr = %.4g) nrm2r = %.4g\n",
638  nn, eigr, eigi, 1. / eigr, nrm2r);
639  std::ostringstream o1;
640  o1 << "eig_" << nn << ".h5m";
641  CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
642  "ELASTIC_MECHANICS", "SPATIAL_POSITION", "EIGEN_VECTOR", COL, D,
643  INSERT_VALUES, SCATTER_REVERSE);
644  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
645  post_proc);
646  CHKERR post_proc.writeFile(o1.str().c_str());
647  }
648 
649  CHKERR KSPDestroy(&solver);
650  CHKERR VecDestroy(&F);
651  CHKERR VecDestroy(&D);
652  CHKERR MatDestroy(&Aij);
653  CHKERR MatDestroy(&Bij);
654  CHKERR EPSDestroy(&eps);
655  }
656  CATCH_ERRORS;
657 
658  SlepcFinalize();
659  MoFEM::Core::Finalize();
660  // PetscFinalize();
661 
662  return 0;
663 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
SIDESET
@ SIDESET
Definition: definitions.h:160
MoFEM::CoreInterface::problem_basic_method_postProcess
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:87
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MyMat_double::D_mu
MatrixDouble D_mu
Definition: stability.cpp:28
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NeumannForcesSurfaceComplexForLazy
NonLinear surface pressure element (obsolete implementation)
Definition: SurfacePressureComplexForLazy.hpp:16
MetaNodalForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI
Implementation of elastic (non-linear) St. Kirchhoff equation.
Definition: NonLinearElasticElement.hpp:79
MyMat_double
Definition: stability.cpp:21
PRESSURESET
@ PRESSURESET
Definition: definitions.h:165
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
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::SnesMethod::snes_x
Vec & snes_x
state vector
Definition: LoopMethods.hpp:121
MoFEM::SnesMethod::snes_ctx
SNESContext snes_ctx
Definition: LoopMethods.hpp:118
MyMat_double::i
FTensor::Index< 'i', 3 > i
Definition: stability.cpp:31
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MyMat::setUserActiveVariables
virtual MoFEMErrorCode setUserActiveVariables(VectorDouble &active_variables)
Add additional independent variables More complex physical models depend on gradient of defamation an...
Definition: stability.cpp:141
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
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
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addForce
MoFEMErrorCode addForce(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:516
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
NodalForce
Force applied to nodes.
Definition: NodalForce.hpp:13
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
MoFEM::SnesMethod::CTX_SNESSETJACOBIAN
@ CTX_SNESSETJACOBIAN
Definition: LoopMethods.hpp:107
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::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::SnesMethod::CTX_SNESSETFUNCTION
@ CTX_SNESSETFUNCTION
Definition: LoopMethods.hpp:107
ROW
@ ROW
Definition: definitions.h:136
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1588
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MetaNodalForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
MyMat
Definition: stability.cpp:111
NODESET
@ NODESET
Definition: definitions.h:159
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
MoFEM::CoreInterface::problem_basic_method_preProcess
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
SurfacePressureComplexForLazy.hpp
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
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
FORCESET
@ FORCESET
Definition: definitions.h:164
MyMat_double::doAotherwiseB
bool doAotherwiseB
Definition: stability.cpp:25
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MyMat::nbActiveVariables0
int nbActiveVariables0
Definition: stability.cpp:113
NonlinearElasticElement::BlockData::E
double E
Definition: HookeElement.hpp:34
Hooke.hpp
Implementation of linear elastic material.
convert.type
type
Definition: convert.py:64
MyMat_double::k
FTensor::Index< 'k', 3 > k
Definition: stability.cpp:33
MyMat_double::calculateP_PiolaKirchhoffI
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
Definition: stability.cpp:35
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
COL
@ COL
Definition: definitions.h:136
help
static char help[]
Definition: stability.cpp:18
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
MyMat_double::MyMat_double
MyMat_double()
Definition: stability.cpp:26
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::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:32
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE
Definition: SurfacePressureComplexForLazy.hpp:39
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::ProblemsManager::partitionFiniteElements
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
Definition: ProblemsManager.cpp:2167
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
NonlinearElasticElement::BlockData::PoissonRatio
double PoissonRatio
Definition: HookeElement.hpp:35
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
FTensor::Index< 'i', 3 >
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
main
int main(int argc, char *argv[])
Definition: stability.cpp:165
NeumannForcesSurfaceComplexForLazy::getLoopSpatialFe
MyTriangleSpatialFE & getLoopSpatialFe()
Definition: SurfacePressureComplexForLazy.hpp:163
MoFEM::Types::VectorBoundedArray
ublas::vector< T, ublas::bounded_array< T, N > > VectorBoundedArray
Definition: Types.hpp:85
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
Range
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
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MyMat_double::sTress
VectorBoundedArray< TYPE, 6 > sTress
Definition: stability.cpp:29
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MyMat::setUserActiveVariables
virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables)
add additional active variables
Definition: stability.cpp:115
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::ProblemsManager::buildProblemOnDistributedMesh
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
Definition: ProblemsManager.cpp:290
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addPressure
MoFEMErrorCode addPressure(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:531
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CoreInterface::modify_problem_add_finite_element
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
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
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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.
MyMat_double::j
FTensor::Index< 'j', 3 > j
Definition: stability.cpp:32
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2339
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
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.
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1683
tol
double tol
Definition: mesh_smoothing.cpp:26
MU
#define MU(E, NU)
Definition: fem_tools.h:23
F
@ F
Definition: free_surface.cpp:394
DirichletSpatialPositionsBc
Set Dirichlet boundary conditions on spatial displacements.
Definition: DirichletBC.hpp:211
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow