v0.12.1
cornea_z_approx.cpp
Go to the documentation of this file.
1 /**
2  * \file cornea_z_approx.cpp
3  * \Cornea cornea_z_approx.cpp
4  *
5  * Using PipelineManager interface calculate the divergence of base functions,
6  * and integral of flux on the boundary. Since the h-div space is used, volume
7  * integral and boundary integral should give the same result.
8  */
9 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
30 #include <BasicFiniteElements.hpp>
31 
32 constexpr int FIELD_DIM = 1;
33 constexpr int SPACE_DIM = 2;
34 
35 template <int DIM> struct ElementsAndOps {};
36 
37 template <> struct ElementsAndOps<2> {
41 };
42 
46 
47 using ThicknessFun = boost::function<double(const double, const double)>;
48 
49 constexpr double Rs = 7.537;
50 constexpr double Rf = 7.755;
51 constexpr double Qs = 0;
52 constexpr double Qf = 0;
53 constexpr double z0 = 2.486;
54 constexpr double alpha = 80 * M_PI / 180;
55 
56 constexpr double thickness_apex = 0.561;
57 constexpr double thickness_limbus = 0.669;
58 constexpr double graft_thickness = 0.1;
60 
61  ApproxFieldFunctionGeneric(ThicknessFun thickness) : zThickness(thickness) {}
62 
63  double operator()(const double x, const double y, const double z) {
64 
65  const double r = sqrt(x * x + y * y);
66  const double theta = atan2(x, y);
67  const double A =
68  (pow(cos(theta - alpha), 2)) / Rs + (pow(sin(theta - alpha), 2)) / Rf;
69  const double B = (Qs + 1) * (pow(cos(theta - alpha), 2)) / (pow(Rs, 2)) +
70  (Qf + 1) * (pow(sin(theta - alpha), 2)) / (pow(Rf, 2));
71  const double calc_z =
72  z0 - (pow(r, 2) * A) / (1 + sqrt(1 - (pow(r, 2) * B)));
73 
74  return calc_z + zThickness(r, theta);
75  }
76 
77 private:
79 };
80 
81 auto thickness_bottom = [](double r, const double theta) {
82  constexpr double R = (Rf + Rs) / 2;
83  const double alpha = 1 - r / R;
84  return -(alpha * thickness_apex + (1 - alpha) * thickness_limbus);
85 };
86 
87 auto thickness_gap = [](double r, const double theta) {
88  return thickness_bottom(r, theta);
89 };
90 
91 auto thickness_graft = [](double r, const double theta) {
92  return thickness_gap(r, theta) + graft_thickness;
93 };
94 
99 
100 struct Cornea {
101 
102  Cornea(MoFEM::Interface &m_field) : mField(m_field) {}
103 
104  MoFEMErrorCode runProblem();
105 
106 private:
108 
109  Range entsTop;
110  Range entsBottom;
111  Range entsGap;
112  Range entsGraft;
113  Range prismsGap;
114  Range prismsGraft;
115 
116  MoFEMErrorCode getOptions();
117  MoFEMErrorCode readMesh();
118  MoFEMErrorCode setupProblem();
119  MoFEMErrorCode setIntegrationRules();
120  MoFEMErrorCode createCommonData();
121  MoFEMErrorCode boundaryCondition();
122  MoFEMErrorCode setOPs();
123  MoFEMErrorCode solveSystem();
124  MoFEMErrorCode outputResults();
125  MoFEMErrorCode checkResults();
126  MoFEMErrorCode setHOGeometry();
127 
128  struct CommonData {
129  boost::shared_ptr<VectorDouble> approxVals;
132  };
133  boost::shared_ptr<CommonData> commonDataPtr;
134 
135  char meshFileName[255]; ///< mesh file name
137  dM; ///< Discrete manager (interface to PETSc/MoFEM functions)
138  boost::shared_ptr<DomainEle> feZRhsTopPtr; ///< FE instance
139  boost::shared_ptr<DomainEle> feZRhsBottomPtr; ///< FE instance
140  boost::shared_ptr<DomainEle> feZRhsGapPtr; ///< FE instance
141  boost::shared_ptr<DomainEle> feZRhsGraftPtr; ///< FE instance
142  boost::shared_ptr<DomainEle> feZLhsPtr; ///< FE instance
143 
144  struct OpError;
145 };
146 
147 struct Cornea::OpError : public DomainEleOp {
148 
149  OpError(boost::shared_ptr<CommonData> &common_data_ptr, ScalarFun z_fun)
150  : DomainEleOp("Z", OPROW), commonDataPtr(common_data_ptr), zFun(z_fun) {}
151  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
153 
154  if (const size_t nb_dofs = data.getIndices().size()) {
155 
156  const int nb_integration_pts = getGaussPts().size2();
157  auto t_w = getFTensor0IntegrationWeight();
158  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
159  auto t_coords = getFTensor1CoordsAtGaussPts();
160 
161  VectorDouble nf(nb_dofs, false);
162  nf.clear();
163 
165  const double volume = getMeasure();
166 
167  auto t_row_base = data.getFTensor0N();
168  double error = 0;
169  for (int gg = 0; gg != nb_integration_pts; ++gg) {
170 
171  const double alpha = t_w * volume;
172  double diff = t_val - zFun(t_coords(0), t_coords(1), t_coords(2));
173  error += alpha * pow(diff, 2);
174 
175  for (size_t r = 0; r != nb_dofs; ++r) {
176  nf[r] += alpha * t_row_base * diff;
177  ++t_row_base;
178  }
179 
180  ++t_w;
181  ++t_val;
182  ++t_coords;
183  }
184 
185  const int index = 0;
186  CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
187  CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
188  }
189 
191  }
192 
193 private:
194  boost::shared_ptr<CommonData> commonDataPtr;
196 };
197 
198 //! [Run programme]
201  CHKERR getOptions();
202  CHKERR readMesh();
203  CHKERR setupProblem();
204  CHKERR setIntegrationRules();
205  CHKERR createCommonData();
206  CHKERR boundaryCondition();
207  CHKERR setOPs();
208  CHKERR solveSystem();
209  CHKERR setHOGeometry();
210  CHKERR checkResults();
211  CHKERR outputResults();
213 }
214 //! [Run programme]
215 
217  PetscBool flg = PETSC_TRUE;
219  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Cornea", "none");
220  CHKERRG(ierr);
221  ierr = PetscOptionsString("-file_name", "file name", "", "mesh.h5m",
222  meshFileName, 255, &flg);
223  CHKERRG(ierr);
224  ierr = PetscOptionsEnd();
225  CHKERRG(ierr);
227 }
228 
229 //! [Read mesh]
232 
233  // This is a case of distributed mesh and algebra. In that case each processor
234  // keep only part of the problem.
235  CHKERR mField.get_moab().load_file(meshFileName, 0, "");
236  CHKERR mField.rebuild_database();
237 
238  auto getBlokset = [&](auto block_name) {
239  const CubitMeshSets *blockset_ptr;
240  CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
241  block_name, &blockset_ptr);
242  return blockset_ptr;
243  };
244 
245  auto getEnts = [&](auto block_name) {
246  auto blockset_ptr = getBlokset(block_name);
247  Range faces;
248  CHKERR blockset_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
249  faces, true);
250  return faces;
251  };
252 
253  entsTop = getEnts("TOP");
254  entsBottom = getEnts("BOTTOM");
255 
256  auto add_two_layers_of_prisms = [&]() {
258  auto interface_ptr = mField.getInterface<PrismsFromSurfaceInterface>();
259  CHKERR interface_ptr->createPrisms(entsBottom, prismsGap);
260  std::array<double, 3> d3{0, 0, 0};
261  std::array<double, 3> d4{0, 0, -1};
262  CHKERR interface_ptr->setThickness(prismsGap, d3.data(), d4.data());
263 
264  interface_ptr->createdVertices.clear();
265  CHKERR interface_ptr->createPrismsFromPrisms(prismsGap, false, prismsGraft);
266  CHKERR interface_ptr->setThickness(prismsGraft, d3.data(), d4.data());
267 
268  auto meshset_mng = mField.getInterface<MeshsetsManager>();
269 
270  CHKERR meshset_mng->addMeshset(BLOCKSET, 101, "GAP");
271  CHKERR meshset_mng->addMeshset(BLOCKSET, 102, "GRAFT");
272  CHKERR meshset_mng->addEntitiesToMeshset(BLOCKSET, 101, prismsGap);
273  CHKERR meshset_mng->addEntitiesToMeshset(BLOCKSET, 102, prismsGraft);
274 
275  Skinner skin(&mField.get_moab());
276  CHKERR skin.find_skin(0, prismsGap, false, entsGap);
277  entsGap = subtract(entsGap, entsBottom).subset_by_type(MBTRI);
278  CHKERR skin.find_skin(0, prismsGraft, false, entsGraft);
279  entsGraft = subtract(entsGraft, entsGap).subset_by_type(MBTRI);
280 
281  CHKERR meshset_mng->addMeshset(BLOCKSET, 103, "SKIN_GAP");
282  CHKERR meshset_mng->addMeshset(BLOCKSET, 104, "SKIN_GRAFT");
283  CHKERR meshset_mng->addEntitiesToMeshset(BLOCKSET, 103, entsGap);
284  CHKERR meshset_mng->addEntitiesToMeshset(BLOCKSET, 104, entsGraft);
285 
287  };
288 
289  CHKERR add_two_layers_of_prisms();
290 
291  Range ents;
292  CHKERR mField.get_moab().get_entities_by_dimension(0, 3, ents, true);
293  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevel(
294  ents, BitRefLevel().set(0), false);
296 }
297 //! [Read mesh]
298 
299 //! [Set up problem]
302 
303  // Declare fields
304  CHKERR mField.add_field("Z", H1, AINSWORTH_LEGENDRE_BASE, 1);
305  CHKERR mField.add_field("HO_GEOMETRY", H1, AINSWORTH_LEGENDRE_BASE, 3);
306 
307  // Set order and add entities to field
308  constexpr int order = 4;
309  Range z_ents = unite(entsTop, entsBottom);
310  z_ents.merge(entsGap);
311  z_ents.merge(entsGraft);
312  CHKERR mField.add_ents_to_field_by_dim(z_ents, 2, "Z");
313  CHKERR mField.set_field_order(0, MBVERTEX, "Z", 1);
314  CHKERR mField.set_field_order(0, MBEDGE, "Z", order);
315  CHKERR mField.set_field_order(0, MBTRI, "Z", order);
316 
317  CHKERR mField.add_ents_to_field_by_type(0, MBTRI, "HO_GEOMETRY");
318  CHKERR mField.set_field_order(0, MBVERTEX, "HO_GEOMETRY", 1);
319  CHKERR mField.set_field_order(0, MBEDGE, "HO_GEOMETRY", order);
320  CHKERR mField.set_field_order(0, MBTRI, "HO_GEOMETRY", order);
321 
322  auto add_ele = [&](auto name, auto ents) {
324  CHKERR mField.add_finite_element(name);
325  // Declare finite ele,elements
326  CHKERR mField.modify_finite_element_add_field_row(name, "Z");
327  CHKERR mField.modify_finite_element_add_field_col(name, "Z");
328  CHKERR mField.modify_finite_element_add_field_data(name, "Z");
329  CHKERR mField.modify_finite_element_add_field_data(name, "HO_GEOMETRY");
330  // Add entities to finite element
331  CHKERR mField.add_ents_to_finite_element_by_dim(ents, 2, name);
333  };
334 
335  CHKERR add_ele("Z_TOP", entsTop);
336  CHKERR add_ele("Z_BOTTOM", entsBottom);
337  CHKERR add_ele("Z_GAP", entsGap);
338  CHKERR add_ele("Z_GRAFT", entsGraft);
339 
340  // Build fields and finite elements
341  CHKERR mField.build_fields();
342  CHKERR mField.build_finite_elements("Z_TOP");
343  CHKERR mField.build_finite_elements("Z_BOTTOM");
344  CHKERR mField.build_finite_elements("Z_GAP");
345  CHKERR mField.build_finite_elements("Z_GRAFT");
346  CHKERR mField.build_adjacencies(BitRefLevel().set(0));
347 
348  dM = createSmartDM(PETSC_COMM_WORLD, "DMMOFEM");
349  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "Z", BitRefLevel().set(0));
350  CHKERR DMSetFromOptions(dM);
351  CHKERR DMMoFEMAddElement(dM, "Z_TOP");
352  CHKERR DMMoFEMAddElement(dM, "Z_BOTTOM");
353  CHKERR DMMoFEMAddElement(dM, "Z_GAP");
354  CHKERR DMMoFEMAddElement(dM, "Z_GRAFT");
355  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
356  CHKERR DMSetUp_MoFEM(dM);
357 
359 }
360 //! [Set up problem]
361 
362 //! [Set integration rule]
365  auto make_fe_ptr = [&](auto &fe_ptr) {
366  fe_ptr = boost::make_shared<DomainEle>(mField);
367  fe_ptr->getRuleHook = [](int, int, int p) -> int { return 2 * p; };
368  };
369  make_fe_ptr(feZRhsTopPtr);
370  make_fe_ptr(feZRhsBottomPtr);
371  make_fe_ptr(feZRhsGapPtr);
372  make_fe_ptr(feZRhsGraftPtr);
373  make_fe_ptr(feZLhsPtr);
375 }
376 //! [Set integration rule]
377 
378 //! [Create common data]
381  commonDataPtr = boost::make_shared<CommonData>();
382  commonDataPtr->resVec = smartCreateDMVector(dM);
383  commonDataPtr->L2Vec = createSmartVectorMPI(
384  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
385  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
387 }
388 //! [Create common data]
389 
390 //! [Boundary condition]
392 //! [Boundary condition]
393 
394 //! [Push operators to pipeline]
397  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
398  auto beta = [](const double, const double, const double) { return 1; };
399  feZLhsPtr->getOpPtrVector().push_back(new OpDomainMass("Z", "Z", beta));
400  feZRhsTopPtr->getOpPtrVector().push_back(new OpDomainSource(
401  "Z", ApproxFieldFunctionGeneric([](double, double) { return 0; })));
402  feZRhsBottomPtr->getOpPtrVector().push_back(
404  feZRhsGapPtr->getOpPtrVector().push_back(
406  feZRhsGraftPtr->getOpPtrVector().push_back(
409 }
410 //! [Push operators to pipeline]
411 
412 //! [Solve]
415 
416  auto solver = MoFEM::createKSP(mField.get_comm());
417  CHKERR KSPSetDM(solver, dM);
418 
419  boost::shared_ptr<FEMethod> null;
420  auto add_fe_instance = [&](auto fe_name, auto fe_lhs, auto fe_rhs) {
422  CHKERR DMMoFEMKSPSetComputeOperators(dM, fe_name, fe_lhs, null, null);
423  CHKERR DMMoFEMKSPSetComputeRHS(dM, fe_name, fe_rhs, null, null);
425  };
426  CHKERR add_fe_instance("Z_TOP", feZLhsPtr, feZRhsTopPtr);
427  CHKERR add_fe_instance("Z_BOTTOM", feZLhsPtr, feZRhsBottomPtr);
428  CHKERR add_fe_instance("Z_GAP", feZLhsPtr, feZRhsGapPtr);
429  CHKERR add_fe_instance("Z_GRAFT", feZLhsPtr, feZRhsGraftPtr);
430  CHKERR KSPSetFromOptions(solver);
431  CHKERR KSPSetUp(solver);
432 
433  auto D = smartCreateDMVector(dM);
434  auto F = smartVectorDuplicate(D);
435 
436  CHKERR KSPSolve(solver, F, D);
437  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
438  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
439  CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_REVERSE);
441 }
442 
443 //! [Solve]
446  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
447  post_proc_fe->generateReferenceElementMesh();
448  post_proc_fe->addFieldValuesPostProc("Z");
449  post_proc_fe->addFieldValuesPostProc("HO_GEOMETRY");
450  CHKERR DMoFEMLoopFiniteElements(dM, "Z_TOP", post_proc_fe);
451  CHKERR DMoFEMLoopFiniteElements(dM, "Z_BOTTOM", post_proc_fe);
452  CHKERR DMoFEMLoopFiniteElements(dM, "Z_GAP", post_proc_fe);
453  CHKERR DMoFEMLoopFiniteElements(dM, "Z_GRAFT", post_proc_fe);
454  CHKERR post_proc_fe->writeFile("out_approx.h5m");
455 
456  dM.reset();
457  CHKERR mField.delete_finite_element("Z_TOP");
458  CHKERR mField.delete_finite_element("Z_BOTTOM");
459  CHKERR mField.delete_finite_element("Z_GAP");
460  CHKERR mField.delete_finite_element("Z_GRAFT");
461  CHKERR mField.delete_field("Z");
462 
463  CHKERR mField.get_moab().write_file("out_ho_mesh.h5m");
464 
466 }
467 //! [Postprocess results]
468 
469 //! [Check results]
472 
473  auto create_check_fe = [&](auto approx_fun) {
474  auto fe = boost::make_shared<DomainEle>(mField);
475  fe->getOpPtrVector().push_back(
476  new OpCalculateScalarFieldValues("Z", commonDataPtr->approxVals));
477  fe->getOpPtrVector().push_back(new OpError(commonDataPtr, approx_fun));
478  return fe;
479  };
480 
481  CHKERR DMoFEMLoopFiniteElements(dM, "Z_TOP",
482  create_check_fe(ApproxFieldFunctionGeneric(
483  [](double, double) { return 0; })));
485  dM, "Z_BOTTOM",
486  create_check_fe(ApproxFieldFunctionGeneric(thickness_bottom)));
487 
488  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
489  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
490  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
491  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
492  double nrm2;
493  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
494  const double *array;
495  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
496  if (mField.get_comm_rank() == 0)
497  PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", sqrt(array[0]),
498  nrm2);
499  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
500 
502 }
503 //! [Check results]
504 
507 
508  EntityMethod set_coords_method;
509  set_coords_method.preProcessHook = []() { return 0; };
510  set_coords_method.operatorHook = [&]() {
512  VectorDouble3 coords(3);
513  EntityHandle v = set_coords_method.entPtr->getEnt();
514  CHKERR mField.get_moab().get_coords(&v, 1, &coords[0]);
515  for (size_t i = 0; i != 3; ++i)
516  set_coords_method.entPtr->getEntFieldData()[i] = coords[i];
518  };
519  set_coords_method.postProcessHook = []() { return 0; };
520 
521  Range verts;
522  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
523  CHKERR mField.loop_entities("HO_GEOMETRY", set_coords_method, &verts);
524 
525  Range ents_top_and_bottom = unite(entsTop, entsBottom);
526  ents_top_and_bottom.merge(entsGap);
527  ents_top_and_bottom.merge(entsGraft);
528  Range sub_ents_top_and_bottom;
529  for (size_t d = 0; d != 2; ++d)
530  CHKERR mField.get_moab().get_adjacencies(ents_top_and_bottom, d, false,
531  sub_ents_top_and_bottom,
532  moab::Interface::UNION);
533  ents_top_and_bottom.merge(sub_ents_top_and_bottom);
534 
535  map<EntityHandle, VectorDouble> z_vec_map;
536  EntityMethod get_z_vec;
537  get_z_vec.preProcessHook = []() { return 0; };
538  get_z_vec.operatorHook = [&]() {
540  EntityHandle v = get_z_vec.entPtr->getEnt();
541  z_vec_map[v] = get_z_vec.entPtr->getEntFieldData();
543  };
544  get_z_vec.postProcessHook = []() { return 0; };
545  CHKERR mField.loop_entities("Z", get_z_vec, &ents_top_and_bottom);
546 
547  EntityMethod set_z_vec;
548  set_z_vec.preProcessHook = []() { return 0; };
549  set_z_vec.operatorHook = [&]() {
551  EntityHandle v = set_z_vec.entPtr->getEnt();
552  int i = 2;
553  for (auto z_dof : z_vec_map.at(v)) {
554  set_z_vec.entPtr->getEntFieldData()[i] = z_dof;
555  i += 3;
556  }
558  };
559  set_z_vec.postProcessHook = []() { return 0; };
560  CHKERR mField.loop_entities("HO_GEOMETRY", set_z_vec, &ents_top_and_bottom);
562 }
563 
564 int main(int argc, char *argv[]) {
565 
566  // Initialisation of MoFEM/PETSc and MOAB data structures
567  const char param_file[] = "param_file.petsc";
568  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
569 
570  try {
571 
572  //! [Register MoFEM discrete manager in PETSc]
573  DMType dm_name = "DMMOFEM";
574  CHKERR DMRegister_MoFEM(dm_name);
575  //! [Register MoFEM discrete manager in PETSc
576 
577  //! [Create MoAB]
578  moab::Core mb_instance; ///< mesh database
579  moab::Interface &moab = mb_instance; ///< mesh database interface
580  //! [Create MoAB]
581 
582  //! [Create MoFEM]
583  MoFEM::Core core(moab); ///< finite element database
584  MoFEM::Interface &m_field = core; ///< finite element database insterface
585  //! [Create MoFEM]
586 
587  //! [Cornea]
588  Cornea ex(m_field);
589  CHKERR ex.runProblem();
590  //! [Cornea]
591  }
592  CATCH_ERRORS;
593 
595 }
static Index< 'p', 3 > p
std::string param_file
int main(int argc, char *argv[])
static char help[]
constexpr double Qs
constexpr double Rs
constexpr double thickness_apex
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
auto thickness_gap
constexpr int SPACE_DIM
constexpr double Rf
constexpr int FIELD_DIM
auto thickness_bottom
constexpr double z0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
DataForcesAndSourcesCore::EntData EntData
auto thickness_graft
constexpr double Qf
constexpr double graft_thickness
constexpr double thickness_limbus
boost::function< double(const double, const double)> ThicknessFun
constexpr double alpha
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ BLOCKSET
Definition: definitions.h:161
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:458
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:478
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMMoFEM.cpp:592
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:541
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMMoFEM.cpp:633
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition: DMMMoFEM.cpp:1200
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:103
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMMoFEM.cpp:393
auto createSmartDM
Creates smart DM object.
auto createKSP
auto createSmartVectorMPI
Create MPI Vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1095
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:152
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1949
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
const double r
rate factor
double A
static constexpr std::array< double, 3 > d4
static constexpr std::array< double, 3 > d3
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ApproxFieldFunctionGeneric(ThicknessFun thickness)
double operator()(const double x, const double y, const double z)
SmartPetscObj< Vec > resVec
SmartPetscObj< Vec > L2Vec
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
OpError(boost::shared_ptr< CommonData > &common_data_ptr, ScalarFun z_fun)
boost::shared_ptr< DomainEle > feZRhsGapPtr
FE instance.
SmartPetscObj< DM > dM
Discrete manager (interface to PETSc/MoFEM functions)
boost::shared_ptr< DomainEle > feZRhsTopPtr
FE instance.
boost::shared_ptr< DomainEle > feZRhsBottomPtr
FE instance.
Range entsGraft
MoFEMErrorCode setIntegrationRules()
[Set up problem]
Cornea(MoFEM::Interface &m_field)
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode setHOGeometry()
[Check results]
MoFEMErrorCode setOPs()
[Boundary condition]
MoFEMErrorCode createCommonData()
[Set integration rule]
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve]
boost::shared_ptr< CommonData > commonDataPtr
MoFEM::Interface & mField
Range entsBottom
MoFEMErrorCode boundaryCondition()
[Create common data]
boost::shared_ptr< DomainEle > feZLhsPtr
FE instance.
MoFEMErrorCode getOptions()
[Run programme]
Range prismsGap
MoFEMErrorCode checkResults()
[Postprocess results]
boost::shared_ptr< DomainEle > feZRhsGraftPtr
FE instance.
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Range prismsGraft
PipelineManager::FaceEle DomainEle
boost::function< MoFEMErrorCode()> preProcessHook
Hook function for pre-processing.
boost::function< MoFEMErrorCode()> operatorHook
Hook function for operator.
boost::function< MoFEMErrorCode()> postProcessHook
Hook function for post-processing.
Managing BitRefLevels.
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
this struct keeps basic methods for moab meshset about material and boundary conditions
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Deprecated interface functions.
Data structure to exchange data between mofem and User Loop Methods on entities.
boost::shared_ptr< FieldEntity > entPtr
Interface for managing meshsets containing materials and boundary conditions.
Get value at integration points for scalar field.
PipelineManager interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Postprocess on face.