v0.13.0
MixTransportElement.hpp
Go to the documentation of this file.
1 /** \file MixTransportElement.hpp
2  * \brief Mix implementation of transport element
3  *
4  * \ingroup mofem_mix_transport_elem
5  *
6  */
7 
8 /*
9  * This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef _MIX_TRANPORT_ELEMENT_HPP_
24 #define _MIX_TRANPORT_ELEMENT_HPP_
25 
26 namespace MixTransport {
27 
28 /** \brief Mix transport problem
29  \ingroup mofem_mix_transport_elem
30 
31  Note to solve this system you need to use direct solver or proper
32  preconditioner for saddle problem.
33 
34  It is based on \cite arnold2006differential \cite arnold2012mixed \cite
35  reviartthomas1996
36  <https://www.researchgate.net/profile/Richard_Falk/publication/226454406_Differential_Complexes_and_Stability_of_Finite_Element_Methods_I._The_de_Rham_Complex/links/02e7e5214f0426ff77000000.pdf>
37 
38  General problem have form,
39  \f[
40  \mathbf{A} \boldsymbol\sigma + \textrm{grad}[u] = \mathbf{0} \; \textrm{on} \;
41  \Omega \\ \textrm{div}[\boldsymbol\sigma] = f \; \textrm{on} \; \Omega \f]
42 
43 */
45 
47 
48  /**
49  * \brief definition of volume element
50 
51  * It is used to calculate volume integrals. On volume element we set-up
52  * operators to calculate components of matrix and vector.
53 
54  */
58  int getRule(int order) { return 2 * order + 1; };
59  };
60 
61  MyVolumeFE feVol; ///< Instance of volume element
62 
63  /** \brief definition of surface element
64 
65  * It is used to calculate surface integrals. On volume element are operators
66  * evaluating natural boundary conditions.
67 
68  */
72  int getRule(int order) { return 2 * order + 1; };
73  };
74 
75  MyTriFE feTri; ///< Instance of surface element
76 
77  /**
78  * \brief construction of this data structure
79  */
81  : mField(m_field), feVol(m_field), feTri(m_field){};
82 
83  /**
84  * \brief destructor
85  */
86  virtual ~MixTransportElement() {}
87 
88  VectorDouble valuesAtGaussPts; ///< values at integration points on element
90  valuesGradientAtGaussPts; ///< gradients at integration points on element
92  divergenceAtGaussPts; ///< divergence at integration points on element
93  MatrixDouble fluxesAtGaussPts; ///< fluxes at integration points on element
94 
95  set<int> bcIndices;
96 
97  /**
98  * \brief get dof indices where essential boundary conditions are applied
99  * @param is indices
100  * @return error code
101  */
104  std::vector<int> ids;
105  ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
106  IS is_local;
107  CHKERR ISCreateGeneral(mField.get_comm(), ids.size(),
108  ids.empty() ? PETSC_NULL : &ids[0],
109  PETSC_COPY_VALUES, &is_local);
110  CHKERR ISAllGather(is_local, is);
111  CHKERR ISDestroy(&is_local);
113  }
114 
115  /**
116  * \brief set source term
117  * @param ent handle to entity on which function is evaluated
118  * @param x coord
119  * @param y coord
120  * @param z coord
121  * @param flux reference to source term set by function
122  * @return error code
123  */
124  virtual MoFEMErrorCode getSource(const EntityHandle ent, const double x,
125  const double y, const double z,
126  double &flux) {
128  flux = 0;
130  }
131 
132  /**
133  * \brief natural (Dirichlet) boundary conditions (set values)
134  * @param ent handle to finite element entity
135  * @param x coord
136  * @param y coord
137  * @param z coord
138  * @param value reference to value set by function
139  * @return error code
140  */
141  virtual MoFEMErrorCode getResistivity(const EntityHandle ent, const double x,
142  const double y, const double z,
143  MatrixDouble3by3 &inv_k) {
145  inv_k.clear();
146  for (int dd = 0; dd < 3; dd++) {
147  inv_k(dd, dd) = 1;
148  }
150  }
151 
152  /**
153  * \brief evaluate natural (Dirichlet) boundary conditions
154  * @param ent entity on which bc is evaluated
155  * @param x coordinate
156  * @param y coordinate
157  * @param z coordinate
158  * @param value vale
159  * @return error code
160  */
161  virtual MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
162  const double x, const double y,
163  const double z, double &value) {
165  value = 0;
167  }
168 
169  /**
170  * \brief essential (Neumann) boundary condition (set fluxes)
171  * @param ent handle to finite element entity
172  * @param x coord
173  * @param y coord
174  * @param z coord
175  * @param flux reference to flux which is set by function
176  * @return [description]
177  */
178  virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
179  const double y, const double z,
180  double &flux) {
182  flux = 0;
184  }
185 
186  /** \brief data for evaluation of het conductivity and heat capacity elements
187  * \ingroup mofem_thermal_elem
188  */
189  struct BlockData {
190  double cOnductivity;
191  double cApacity;
192  Range tEts; ///< constatins elements in block set
193  };
194  std::map<int, BlockData>
195  setOfBlocks; ///< maps block set id with appropriate BlockData
196 
197  /**
198  * \brief Add fields to database
199  * @param values name of the fields
200  * @param fluxes name of filed for fluxes
201  * @param order order of approximation
202  * @return error code
203  */
204  MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
205  const int order) {
207  // Fields
210 
211  // meshset consisting all entities in mesh
212  EntityHandle root_set = mField.get_moab().get_root_set();
213  // add entities to field
214  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
215  CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
216  CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
217  CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
218  CHKERR mField.set_field_order(root_set, MBTET, values, order);
220  }
221 
222  /// \brief add finite elements
224  const std::string &fluxes_name, const std::string &values_name) {
226 
227  // Set up volume element operators. Operators are used to calculate
228  // components of stiffness matrix & right hand side, in essence are used to
229  // do volume integrals over tetrahedral in this case.
230 
231  // Define element "MIX". Note that this element will work with fluxes_name
232  // and values_name. This reflect bilinear form for the problem
240 
241  // Define finite element to integrate over skeleton, we need that to
242  // evaluate error
243  CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
245  fluxes_name);
247  fluxes_name);
249  fluxes_name);
250 
251  // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
252  // BLOCKSETS and add them to "MIX" finite element. In addition get data form
253  // that meshset and set conductivity which is used to calculate fluxes from
254  // gradients of concentration or gradient of temperature, depending how you
255  // interpret variables.
257  mField, BLOCKSET | MAT_THERMALSET, it)) {
258 
259  Mat_Thermal temp_data;
260  CHKERR it->getAttributeDataStructure(temp_data);
261  setOfBlocks[it->getMeshsetId()].cOnductivity =
262  temp_data.data.Conductivity;
263  setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
264  CHKERR mField.get_moab().get_entities_by_type(
265  it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
267  setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
268 
269  Range skeleton;
270  CHKERR mField.get_moab().get_adjacencies(
271  setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
272  moab::Interface::UNION);
274  "MIX_SKELETON");
275  }
276 
277  // Define element to integrate natural boundary conditions, i.e. set values.
278  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
280  fluxes_name);
282  fluxes_name);
284  fluxes_name);
286  values_name);
287 
288  // Define element to apply essential boundary conditions.
289  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
291  fluxes_name);
293  fluxes_name);
295  fluxes_name);
297  values_name);
299  }
300 
301  /**
302  * \brief Build problem
303  * @param ref_level mesh refinement on which mesh problem you like to built.
304  * @return error code
305  */
308  // build field
310  // get tetrahedrons which has been build previously and now in so called
311  // garbage bit level
312  Range done_tets;
313  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
314  BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
315  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
316  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
317  done_tets);
318  // get tetrahedrons which belong to problem bit level
319  Range ref_tets;
320  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
321  ref_level, BitRefLevel().set(), MBTET, ref_tets);
322  ref_tets = subtract(ref_tets, done_tets);
323  CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
324  // get triangles which has been build previously and now in so called
325  // garbage bit level
326  Range done_faces;
327  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
328  BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
329  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
330  BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
331  done_faces);
332  // get triangles which belong to problem bit level
333  Range ref_faces;
334  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
335  ref_level, BitRefLevel().set(), MBTRI, ref_faces);
336  ref_faces = subtract(ref_faces, done_faces);
337  // build finite elements structures
338  CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
339  CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
340  CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
341  // Build adjacencies of degrees of freedom and elements
342  CHKERR mField.build_adjacencies(ref_level);
343  // Define problem
345  // set refinement level for problem
347  // Add element to problem
349  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
350  // Boundary conditions
351  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCFLUX");
352  CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
353  // build problem
354 
355  ProblemsManager *prb_mng_ptr;
356  CHKERR mField.getInterface(prb_mng_ptr);
357  CHKERR prb_mng_ptr->buildProblem("MIX", true);
358  // mesh partitioning
359  // partition
360  CHKERR prb_mng_ptr->partitionProblem("MIX");
361  CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
362  // what are ghost nodes, see Petsc Manual
363  CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
365  }
366 
367  struct OpPostProc
370  std::vector<EntityHandle> &mapGaussPts;
371  OpPostProc(moab::Interface &post_proc_mesh,
372  std::vector<EntityHandle> &map_gauss_pts)
374  "VALUES", UserDataOperator::OPCOL),
375  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
376  MoFEMErrorCode doWork(int side, EntityType type,
379  if (type != MBTET)
381  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
382  Tag th_error_flux;
383  CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_FLUX",
384  th_error_flux);
385  double *error_flux_ptr;
386  CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
387  th_error_flux, &fe_ent, 1, (const void **)&error_flux_ptr);
388 
389  Tag th_error_div;
390  CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_DIV",
391  th_error_div);
392  double *error_div_ptr;
393  CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
394  th_error_div, &fe_ent, 1, (const void **)&error_div_ptr);
395 
396  Tag th_error_jump;
397  CHKERR getVolumeFE()->mField.get_moab().tag_get_handle("ERROR_JUMP",
398  th_error_jump);
399  double *error_jump_ptr;
400  CHKERR getVolumeFE()->mField.get_moab().tag_get_by_ptr(
401  th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
402 
403  {
404  double def_val = 0;
405  Tag th_error_flux;
406  CHKERR postProcMesh.tag_get_handle(
407  "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
408  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
409  for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
410  vit != mapGaussPts.end(); vit++) {
411  CHKERR postProcMesh.tag_set_data(th_error_flux, &*vit, 1,
412  error_flux_ptr);
413  }
414 
415  Tag th_error_div;
416  CHKERR postProcMesh.tag_get_handle(
417  "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
418  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
419  for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
420  vit != mapGaussPts.end(); vit++) {
421  CHKERR postProcMesh.tag_set_data(th_error_div, &*vit, 1,
422  error_div_ptr);
423  }
424 
425  Tag th_error_jump;
426  CHKERR postProcMesh.tag_get_handle(
427  "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
428  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
429  for (vector<EntityHandle>::iterator vit = mapGaussPts.begin();
430  vit != mapGaussPts.end(); vit++) {
431  CHKERR postProcMesh.tag_set_data(th_error_jump, &*vit, 1,
432  error_jump_ptr);
433  }
434  }
436  }
437  };
438 
439  /**
440  * \brief Post process results
441  * @return error code
442  */
443  MoFEMErrorCode postProc(const string out_file) {
447  CHKERR post_proc.addFieldValuesPostProc("VALUES");
448  CHKERR post_proc.addFieldValuesGradientPostProc("VALUES");
449  CHKERR post_proc.addFieldValuesPostProc("FLUXES");
450  // CHKERR post_proc.addHdivFunctionsPostProc("FLUXES");
451  post_proc.getOpPtrVector().push_back(
452  new OpPostProc(post_proc.postProcMesh, post_proc.mapGaussPts));
453  CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
454  CHKERR post_proc.writeFile(out_file.c_str());
456  }
457 
458  Vec D, D0, F;
459  Mat Aij;
460 
461  /// \brief create matrices
464  CHKERR mField.getInterface<MatrixManager>()
465  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
466  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
467  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
468  CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
470  }
471 
472  /**
473  * \brief solve problem
474  * @return error code
475  */
478  CHKERR MatZeroEntries(Aij);
479  CHKERR VecZeroEntries(F);
480  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
481  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
482  CHKERR VecZeroEntries(D0);
483  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
484  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
485  CHKERR VecZeroEntries(D);
486  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
487  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
488 
489  CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
490  "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
491 
492  // Calculate essential boundary conditions
493 
494  // clear essential bc indices, it could have dofs from other mesh refinement
495  bcIndices.clear();
496  // clear operator, just in case if some other operators are left on this
497  // element
498  feTri.getOpPtrVector().clear();
499  feTri.getOpPtrVector().push_back(
500  new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
501  // set operator to calculate essential boundary conditions
502  feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this, "FLUXES"));
503  CHKERR mField.loop_finite_elements("MIX", "MIX_BCFLUX", feTri);
504  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
505  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
506  CHKERR VecAssemblyBegin(D0);
507  CHKERR VecAssemblyEnd(D0);
508 
509  // set operators to calculate matrix and right hand side vectors
510  feVol.getOpPtrVector().clear();
511  feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
512  feVol.getOpPtrVector().push_back(
513  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
514  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
515  feVol.getOpPtrVector().push_back(
516  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
517  feVol.getOpPtrVector().push_back(
518  new OpTauDotSigma_HdivHdiv(*this, "FLUXES", Aij, F));
519  feVol.getOpPtrVector().push_back(
520  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", Aij, F));
521  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
522 
523  // calculate right hand side for natural boundary conditions
524  feTri.getOpPtrVector().clear();
525  feTri.getOpPtrVector().push_back(
526  new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
527  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
528  CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
529 
530  // assemble matrices
531  CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
532  CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
533  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
534  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
535  CHKERR VecAssemblyBegin(F);
536  CHKERR VecAssemblyEnd(F);
537 
538  {
539  double nrm2_F;
540  CHKERR VecNorm(F, NORM_2, &nrm2_F);
541  PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
542  }
543 
544  // CHKERR MatMultAdd(Aij,D0,F,F); ;
545 
546  // for ksp solver vector is moved into rhs side
547  // for snes it is left ond the left
548  CHKERR VecScale(F, -1);
549 
550  IS essential_bc_ids;
551  CHKERR getDirichletBCIndices(&essential_bc_ids);
552  CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0, F);
553  CHKERR ISDestroy(&essential_bc_ids);
554 
555  // {
556  // double norm;
557  // CHKERR MatNorm(Aij,NORM_FROBENIUS,&norm); ;
558  // PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
559  // }
560 
561  {
562  double nrm2_F;
563  CHKERR VecNorm(F, NORM_2, &nrm2_F);
564  PetscPrintf(PETSC_COMM_WORLD, "With BC nrm2_F = %6.4e\n", nrm2_F);
565  }
566 
567  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
568  // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
569  // std::string wait;
570  // std::cin >> wait;
571 
572  // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
573  // std::cin >> wait;
574 
575  // Solve
576  KSP solver;
577  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
578  CHKERR KSPSetOperators(solver, Aij, Aij);
579  CHKERR KSPSetFromOptions(solver);
580  CHKERR KSPSetUp(solver);
581  CHKERR KSPSolve(solver, F, D);
582  CHKERR KSPDestroy(&solver);
583 
584  {
585  double nrm2_D;
586  CHKERR VecNorm(D, NORM_2, &nrm2_D);
587  ;
588  PetscPrintf(PETSC_COMM_WORLD, "nrm2_D = %6.4e\n", nrm2_D);
589  }
590  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
591  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
592 
593  // copy data form vector on mesh
594  CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
595  "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
596 
598  }
599 
600  /// \brief calculate residual
603  CHKERR VecZeroEntries(F);
604  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
605  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
606  CHKERR VecAssemblyBegin(F);
607  CHKERR VecAssemblyEnd(F);
608  // calculate residuals
609  feVol.getOpPtrVector().clear();
610  feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
611  feVol.getOpPtrVector().push_back(
612  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
613  feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
614  feVol.getOpPtrVector().push_back(
615  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
616  feVol.getOpPtrVector().push_back(
617  new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
618  feVol.getOpPtrVector().push_back(
619  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
620  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
621  feTri.getOpPtrVector().clear();
622  feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
623  CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
624  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
625  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
626  CHKERR VecAssemblyBegin(F);
627  CHKERR VecAssemblyEnd(F);
628  // CHKERR VecAXPY(F,-1.,D0);
629  // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
630  {
631  std::vector<int> ids;
632  ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
633  std::vector<double> vals(ids.size(), 0);
634  CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
635  INSERT_VALUES);
636  CHKERR VecAssemblyBegin(F);
637  CHKERR VecAssemblyEnd(F);
638  }
639  {
640  double nrm2_F;
641  CHKERR VecNorm(F, NORM_2, &nrm2_F);
642  PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
643  const double eps = 1e-8;
644  if (nrm2_F > eps) {
645  // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
646  // residual");
647  }
648  }
650  }
651 
652  /** \brief Calculate error on elements
653 
654  \todo this functions runs serial, in future should be parallel and work on
655  distributed meshes.
656 
657  */
660  errorMap.clear();
661  sumErrorFlux = 0;
662  sumErrorDiv = 0;
663  sumErrorJump = 0;
664  feTri.getOpPtrVector().clear();
665  feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
666  CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
668  feVol.getOpPtrVector().clear();
669  feVol.getOpPtrVector().push_back(
670  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
671  feVol.getOpPtrVector().push_back(
672  new OpValuesGradientAtGaussPts(*this, "VALUES"));
673  feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
674  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
676  const Problem *problem_ptr;
677  CHKERR mField.get_problem("MIX", &problem_ptr);
678  PetscPrintf(mField.get_comm(),
679  "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
680  "jump^2 = %6.4e error tot^2 = %6.4e\n",
681  problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
684  }
685 
686  /// \brief destroy matrices
689  CHKERR MatDestroy(&Aij);
690  ;
691  CHKERR VecDestroy(&D);
692  ;
693  CHKERR VecDestroy(&D0);
694  ;
695  CHKERR VecDestroy(&F);
696  ;
698  }
699 
700  /**
701  \brief Assemble \f$\int_\mathcal{T} \mathbf{A} \boldsymbol\sigma \cdot
702  \boldsymbol\tau \textrm{d}\mathcal{T}\f$
703 
704  \ingroup mofem_mix_transport_elem
705  */
708 
710  Mat Aij;
712 
714  const std::string flux_name, Mat aij, Vec f)
716  flux_name, flux_name,
718  cTx(ctx), Aij(aij), F(f) {
719  sYmm = true;
720  }
721 
725 
726  /**
727  * \brief Assemble matrix
728  * @param row_side local index of row entity on element
729  * @param col_side local index of col entity on element
730  * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
731  * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
732  * @param row_data data for row
733  * @param col_data data for col
734  * @return error code
735  */
736  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
737  EntityType col_type,
738  EntitiesFieldData::EntData &row_data,
739  EntitiesFieldData::EntData &col_data) {
741  if (Aij == PETSC_NULL)
743  if (row_data.getIndices().size() == 0)
745  if (col_data.getIndices().size() == 0)
747  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
748  int nb_row = row_data.getIndices().size();
749  int nb_col = col_data.getIndices().size();
750  NN.resize(nb_row, nb_col, false);
751  NN.clear();
754  invK.resize(3, 3, false);
755  // get access to resistivity data by tensor rank 2
757  &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
758  &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
759  // Get base functions
760  auto t_n_hdiv_row = row_data.getFTensor1N<3>();
762  int nb_gauss_pts = row_data.getN().size1();
763  for (int gg = 0; gg != nb_gauss_pts; gg++) {
764  // get integration weight and multiply by element volume
765  double w = getGaussPts()(3, gg) * getVolume();
766  const double x = getCoordsAtGaussPts()(gg, 0);
767  const double y = getCoordsAtGaussPts()(gg, 1);
768  const double z = getCoordsAtGaussPts()(gg, 2);
769  // calculate receptivity (invers of conductivity)
770  CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
771  for (int kk = 0; kk != nb_row; kk++) {
773  &col_data.getVectorN<3>(gg)(0, HVEC0),
774  &col_data.getVectorN<3>(gg)(0, HVEC1),
775  &col_data.getVectorN<3>(gg)(0, HVEC2), 3);
776  t_row(j) = w * t_n_hdiv_row(i) * t_inv_k(i, j);
777  for (int ll = 0; ll != nb_col; ll++) {
778  NN(kk, ll) += t_row(j) * t_n_hdiv_col(j);
779  ++t_n_hdiv_col;
780  }
781  ++t_n_hdiv_row;
782  }
783  }
784  Mat a = (Aij != PETSC_NULL) ? Aij : getFEMethod()->ts_B;
785  CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
786  &col_data.getIndices()[0], &NN(0, 0), ADD_VALUES);
787  // matrix is symmetric, assemble other part
788  if (row_side != col_side || row_type != col_type) {
789  transNN.resize(nb_col, nb_row);
790  noalias(transNN) = trans(NN);
791  CHKERR MatSetValues(a, nb_col, &col_data.getIndices()[0], nb_row,
792  &row_data.getIndices()[0], &transNN(0, 0),
793  ADD_VALUES);
794  }
795 
797  }
798 
799  /**
800  * \brief Assemble matrix
801  * @param side local index of row entity on element
802  * @param type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
803  * @param data data for row
804  * @return error code
805  */
806  MoFEMErrorCode doWork(int side, EntityType type,
809  if (F == PETSC_NULL)
811  int nb_row = data.getIndices().size();
812  if (nb_row == 0)
814 
815  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
816  Nf.resize(nb_row);
817  Nf.clear();
820  invK.resize(3, 3, false);
821  Nf.resize(nb_row);
822  Nf.clear();
823  // get access to resistivity data by tensor rank 2
825  &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
826  &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
827  // get base functions
828  auto t_n_hdiv = data.getFTensor1N<3>();
829  auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
830  int nb_gauss_pts = data.getN().size1();
831  for (int gg = 0; gg < nb_gauss_pts; gg++) {
832  double w = getGaussPts()(3, gg) * getVolume();
833  const double x = getCoordsAtGaussPts()(gg, 0);
834  const double y = getCoordsAtGaussPts()(gg, 1);
835  const double z = getCoordsAtGaussPts()(gg, 2);
836  CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
837  for (int ll = 0; ll != nb_row; ll++) {
838  Nf[ll] += w * t_n_hdiv(i) * t_inv_k(i, j) * t_flux(j);
839  ++t_n_hdiv;
840  }
841  ++t_flux;
842  }
843  CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
844 
846  }
847  };
848 
849  /** \brief Assemble \f$ \int_\mathcal{T} u \textrm{div}[\boldsymbol\tau]
850  * \textrm{d}\mathcal{T} \f$
851  */
854 
857 
858  OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row,
859  string val_name_col, Vec f)
861  flux_name_row, val_name_col, UserDataOperator::OPROW),
862  cTx(ctx), F(f) {
863  // this operator is not symmetric setting this variable makes element
864  // operator to loop over element entities (sub-simplices) without
865  // assumption that off-diagonal matrices are symmetric.
866  sYmm = false;
867  }
868 
870 
871  MoFEMErrorCode doWork(int side, EntityType type,
874  if (data.getFieldData().size() == 0)
876  int nb_row = data.getIndices().size();
877  Nf.resize(nb_row);
878  Nf.clear();
879  divVec.resize(data.getN().size2() / 3, 0);
880  if (divVec.size() != data.getIndices().size()) {
881  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
882  "data inconsistency");
883  }
884  int nb_gauss_pts = data.getN().size1();
885 
887  auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
888 
889  int gg = 0;
890  for (; gg < nb_gauss_pts; gg++) {
891  double w = getGaussPts()(3, gg) * getVolume();
892  for(auto &v : divVec) {
893  v = t_base_diff_hdiv(i, i);
894  ++t_base_diff_hdiv;
895  }
896  noalias(Nf) -= w * divVec * cTx.valuesAtGaussPts[gg];
897  }
898  CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
899  ;
900 
902  }
903  };
904 
905  /** \brief \f$ \int_\mathcal{T} \textrm{div}[\boldsymbol\sigma] v
906  * \textrm{d}\mathcal{T} \f$
907  */
910 
912  Mat Aij;
914 
915  /**
916  * \brief Constructor
917  */
918  OpVDivSigma_L2Hdiv(MixTransportElement &ctx, const std::string val_name_row,
919  string flux_name_col, Mat aij, Vec f)
921  val_name_row, flux_name_col,
923  cTx(ctx), Aij(aij), F(f) {
924 
925  // this operator is not symmetric setting this variable makes element
926  // operator to loop over element entities without
927  // assumption that off-diagonal matrices are symmetric.
928  sYmm = false;
929  }
930 
933 
934  /**
935  * \brief Do calculations
936  * @param row_side local index of entity on row
937  * @param col_side local index of entity on column
938  * @param row_type type of row entity
939  * @param col_type type of col entity
940  * @param row_data row data structure carrying information about base
941  * functions, DOFs indices, etc.
942  * @param col_data column data structure carrying information about base
943  * functions, DOFs indices, etc.
944  * @return error code
945  */
946  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
947  EntityType col_type,
948  EntitiesFieldData::EntData &row_data,
949  EntitiesFieldData::EntData &col_data) {
951  if (Aij == PETSC_NULL)
953  if (row_data.getFieldData().size() == 0)
955  if (col_data.getFieldData().size() == 0)
957  int nb_row = row_data.getFieldData().size();
958  int nb_col = col_data.getFieldData().size();
959  NN.resize(nb_row, nb_col);
960  NN.clear();
961  divVec.resize(nb_col, false);
962 
964  auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
965 
966  int nb_gauss_pts = row_data.getN().size1();
967  for (int gg = 0; gg < nb_gauss_pts; gg++) {
968  double w = getGaussPts()(3, gg) * getVolume();
969  for (auto &v : divVec) {
970  v = t_base_diff_hdiv(i, i);
971  ++t_base_diff_hdiv;
972  }
973  noalias(NN) += w * outer_prod(row_data.getN(gg), divVec);
974  }
975  CHKERR MatSetValues(Aij, nb_row, &row_data.getIndices()[0], nb_col,
976  &col_data.getIndices()[0], &NN(0, 0), ADD_VALUES);
977  transNN.resize(nb_col, nb_row);
978  ublas::noalias(transNN) = -trans(NN);
979  CHKERR MatSetValues(Aij, nb_col, &col_data.getIndices()[0], nb_row,
980  &row_data.getIndices()[0], &transNN(0, 0),
981  ADD_VALUES);
983  }
984 
985  MoFEMErrorCode doWork(int side, EntityType type,
988  if (data.getIndices().size() == 0)
990  if (data.getIndices().size() != data.getN().size2()) {
991  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
992  "data inconsistency");
993  }
994  int nb_row = data.getIndices().size();
995  Nf.resize(nb_row);
996  Nf.clear();
997  int nb_gauss_pts = data.getN().size1();
998  for (int gg = 0; gg < nb_gauss_pts; gg++) {
999  double w = getGaussPts()(3, gg) * getVolume();
1000  noalias(Nf) += w * data.getN(gg) * cTx.divergenceAtGaussPts[gg];
1001  }
1002  CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1003  ;
1005  }
1006  };
1007 
1008  /** \brief Calculate source therms, i.e. \f$\int_\mathcal{T} f v
1009  * \textrm{d}\mathcal{T}\f$
1010  */
1011  struct OpL2Source
1015 
1016  OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
1018  val_name, UserDataOperator::OPROW),
1019  cTx(ctx), F(f) {}
1020 
1022  MoFEMErrorCode doWork(int side, EntityType type,
1025  if (data.getFieldData().size() == 0)
1027  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1028  int nb_row = data.getFieldData().size();
1029  Nf.resize(nb_row, false);
1030  Nf.clear();
1031  int nb_gauss_pts = data.getN().size1();
1032  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1033  double w = getGaussPts()(3, gg) * getVolume();
1034  const double x = getCoordsAtGaussPts()(gg, 0);
1035  const double y = getCoordsAtGaussPts()(gg, 1);
1036  const double z = getCoordsAtGaussPts()(gg, 2);
1037  double flux = 0;
1038  CHKERR cTx.getSource(fe_ent, x, y, z, flux);
1039  noalias(Nf) += w * data.getN(gg) * flux;
1040  }
1041  CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1042  ;
1044  }
1045  };
1046 
1047  /**
1048  * \brief calculate \f$ \int_\mathcal{S} {\boldsymbol\tau} \cdot \mathbf{n}u
1049  \textrm{d}\mathcal{S} \f$
1050 
1051  * This terms comes from differentiation by parts. Note that in this Dirichlet
1052  * boundary conditions are natural.
1053 
1054  */
1057 
1060 
1061  /**
1062  * \brief Constructor
1063  */
1064  OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name,
1065  Vec f)
1067  fluxes_name, UserDataOperator::OPROW),
1068  cTx(ctx), F(f) {}
1069 
1071 
1072  /**
1073  * \brief Integrate boundary condition
1074  * @param side local index of entity
1075  * @param type type of entity
1076  * @param data data on entity
1077  * @return error code
1078  */
1079  MoFEMErrorCode doWork(int side, EntityType type,
1082  if (data.getFieldData().size() == 0)
1084  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1085  nF.resize(data.getIndices().size());
1086  nF.clear();
1087  int nb_gauss_pts = data.getN().size1();
1088  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1089  const double x = getCoordsAtGaussPts()(gg, 0);
1090  const double y = getCoordsAtGaussPts()(gg, 1);
1091  const double z = getCoordsAtGaussPts()(gg, 2);
1092 
1093  double value;
1094  CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
1095  ;
1096  double w = getGaussPts()(2, gg) * 0.5;
1097  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1098  noalias(nF) +=
1099  w * prod(data.getVectorN<3>(gg), getNormalsAtGaussPts(gg)) *
1100  value;
1101  } else {
1102  noalias(nF) += w * prod(data.getVectorN<3>(gg), getNormal()) * value;
1103  }
1104  }
1105  Vec f = (F != PETSC_NULL) ? F : getFEMethod()->ts_F;
1106  CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
1107  &nF[0], ADD_VALUES);
1108 
1110  }
1111  };
1112 
1113  /**
1114  * \brief Evaluate boundary conditions on fluxes.
1115  *
1116  * Note that Neumann boundary conditions here are essential. So it is opposite
1117  * what you find in displacement finite element method.
1118  *
1119 
1120  * Here we have to solve for degrees of freedom on boundary such base
1121  functions
1122  * approximate flux.
1123 
1124  *
1125  */
1129  OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
1131  flux_name, UserDataOperator::OPROW),
1132  cTx(ctx) {}
1133 
1137 
1138  MoFEMErrorCode doWork(int side, EntityType type,
1141  if (data.getFieldData().size() == 0)
1143  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1144  int nb_dofs = data.getFieldData().size();
1145  int nb_gauss_pts = data.getN().size1();
1146  if (3 * nb_dofs != static_cast<int>(data.getN().size2())) {
1147  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1148  "wrong number of dofs");
1149  }
1150  NN.resize(nb_dofs, nb_dofs);
1151  Nf.resize(nb_dofs);
1152  NN.clear();
1153  Nf.clear();
1154 
1155  // Get normal vector. Note that when higher order geometry is set, then
1156  // face element could be curved, i.e. normal can be different at each
1157  // integration point.
1158  double *normal_ptr;
1159  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1160  // HO geometry
1161  normal_ptr = &getNormalsAtGaussPts(0)[0];
1162  } else {
1163  // Linear geometry, i.e. constant normal on face
1164  normal_ptr = &getNormal()[0];
1165  }
1166  // set tensor from pointer
1167  FTensor::Tensor1<const double *, 3> t_normal(normal_ptr, &normal_ptr[1],
1168  &normal_ptr[2], 3);
1169  // get base functions
1170  auto t_n_hdiv_row = data.getFTensor1N<3>();
1171 
1172  double nrm2 = 0;
1173 
1174  // loop over integration points
1175  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1176 
1177  // get integration point coordinates
1178  const double x = getCoordsAtGaussPts()(gg, 0);
1179  const double y = getCoordsAtGaussPts()(gg, 1);
1180  const double z = getCoordsAtGaussPts()(gg, 2);
1181 
1182  // get flux on fece for given element handle and coordinates
1183  double flux;
1184  CHKERR cTx.getBcOnFluxes(fe_ent, x, y, z, flux);
1185  ;
1186  // get weight for integration rule
1187  double w = getGaussPts()(2, gg);
1188  if (gg == 0) {
1189  nrm2 = sqrt(t_normal(i) * t_normal(i));
1190  }
1191 
1192  // set tensor of rank 0 to matrix NN elements
1193  // loop over base functions on rows and columns
1194  for (int ll = 0; ll != nb_dofs; ll++) {
1195  // get column on shape functions
1197  &data.getVectorN<3>(gg)(0, HVEC0),
1198  &data.getVectorN<3>(gg)(0, HVEC1),
1199  &data.getVectorN<3>(gg)(0, HVEC2), 3);
1200  for (int kk = 0; kk <= ll; kk++) {
1201  NN(ll, kk) += w * t_n_hdiv_row(i) * t_n_hdiv_col(i);
1202  ++t_n_hdiv_col;
1203  }
1204  // right hand side
1205  Nf[ll] += w * t_n_hdiv_row(i) * t_normal(i) * flux / nrm2;
1206  ++t_n_hdiv_row;
1207  }
1208 
1209  // If HO geometry increment t_normal to next integration point
1210  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1211  ++t_normal;
1212  nrm2 = sqrt(t_normal(i) * t_normal(i));
1213  }
1214  }
1215  // get global dofs indices on element
1216  cTx.bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1217  // factor matrix
1219  // solve local problem
1220  cholesky_solve(NN, Nf, ublas::lower());
1221 
1222  // set solution to vector
1223  CHKERR VecSetOption(cTx.D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1224  CHKERR VecSetValues(cTx.D0, nb_dofs, &*data.getIndices().begin(),
1225  &*Nf.begin(), INSERT_VALUES);
1226  for (int dd = 0; dd != nb_dofs; ++dd)
1227  data.getFieldDofs()[dd]->getFieldData() = Nf[dd];
1228 
1230  }
1231  };
1232 
1233  /**
1234  * \brief Calculate values at integration points
1235  */
1238 
1240 
1242  const std::string val_name = "VALUES")
1244  val_name, UserDataOperator::OPROW),
1245  cTx(ctx) {}
1246 
1247  virtual ~OpValuesAtGaussPts() {}
1248 
1249  MoFEMErrorCode doWork(int side, EntityType type,
1252  if (data.getFieldData().size() == 0)
1254 
1255  int nb_gauss_pts = data.getN().size1();
1256  cTx.valuesAtGaussPts.resize(nb_gauss_pts);
1257  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1258  cTx.valuesAtGaussPts[gg] =
1259  inner_prod(trans(data.getN(gg)), data.getFieldData());
1260  }
1261 
1263  }
1264  };
1265 
1266  /**
1267  * \brief Calculate gradients of values at integration points
1268  */
1271 
1273 
1275  const std::string val_name = "VALUES")
1277  val_name, UserDataOperator::OPROW),
1278  cTx(ctx) {}
1280 
1281  MoFEMErrorCode doWork(int side, EntityType type,
1284  if (data.getFieldData().size() == 0)
1286  int nb_gauss_pts = data.getDiffN().size1();
1287  cTx.valuesGradientAtGaussPts.resize(3, nb_gauss_pts);
1288  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1289  ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1291  noalias(values_grad_at_gauss_pts) =
1292  prod(trans(data.getDiffN(gg)), data.getFieldData());
1293  }
1295  }
1296  };
1297 
1298  /**
1299  * \brief calculate flux at integration point
1300  */
1303 
1305 
1307  const std::string field_name)
1309  field_name, UserDataOperator::OPROW),
1310  cTx(ctx) {}
1311 
1313  MoFEMErrorCode doWork(int side, EntityType type,
1316 
1317  if (data.getFieldData().size() == 0)
1319  int nb_gauss_pts = data.getDiffN().size1();
1320  int nb_dofs = data.getFieldData().size();
1321  cTx.fluxesAtGaussPts.resize(3, nb_gauss_pts);
1322  cTx.divergenceAtGaussPts.resize(nb_gauss_pts);
1323  if (type == MBTRI && side == 0) {
1324  cTx.divergenceAtGaussPts.clear();
1325  cTx.fluxesAtGaussPts.clear();
1326  }
1327  divVec.resize(nb_dofs);
1328 
1330  auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
1331 
1332  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1333  for(auto &v : divVec) {
1334  v = t_base_diff_hdiv(i, i);
1335  ++t_base_diff_hdiv;
1336  }
1337  cTx.divergenceAtGaussPts[gg] += inner_prod(divVec, data.getFieldData());
1338  ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1339  cTx.fluxesAtGaussPts, gg);
1340  flux_at_gauss_pt +=
1341  prod(trans(data.getVectorN<3>(gg)), data.getFieldData());
1342  }
1344  }
1345  };
1346 
1347  map<double, EntityHandle> errorMap;
1349  double sumErrorDiv;
1351 
1352  /** \brief calculate error evaluator
1353  */
1354  struct OpError
1356 
1358 
1359  OpError(MixTransportElement &ctx, const std::string field_name)
1361  field_name, UserDataOperator::OPROW),
1362  cTx(ctx) {}
1363  virtual ~OpError() {}
1364 
1367 
1368  MoFEMErrorCode doWork(int side, EntityType type,
1371  if (type != MBTET)
1373  invK.resize(3, 3, false);
1374  int nb_gauss_pts = data.getN().size1();
1375  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1376  double def_val = 0;
1377  Tag th_error_flux, th_error_div;
1378  CHKERR cTx.mField.get_moab().tag_get_handle(
1379  "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
1380  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1381  double *error_flux_ptr;
1382  CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1383  th_error_flux, &fe_ent, 1, (const void **)&error_flux_ptr);
1384 
1385  CHKERR cTx.mField.get_moab().tag_get_handle(
1386  "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
1387  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1388  double *error_div_ptr;
1389  CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1390  th_error_div, &fe_ent, 1, (const void **)&error_div_ptr);
1391 
1392  Tag th_error_jump;
1393  CHKERR cTx.mField.get_moab().tag_get_handle(
1394  "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1395  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1396  double *error_jump_ptr;
1397  CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1398  th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
1399  *error_jump_ptr = 0;
1400 
1401  /// characteristic size of the element
1402  const double h = pow(getVolume() * 12 / std::sqrt(2), (double)1 / 3);
1403 
1404  for (int ff = 0; ff != 4; ff++) {
1405  EntityHandle face;
1406  CHKERR cTx.mField.get_moab().side_element(fe_ent, 2, ff, face);
1407  double *error_face_jump_ptr;
1408  CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1409  th_error_jump, &face, 1, (const void **)&error_face_jump_ptr);
1410  *error_face_jump_ptr = (1 / sqrt(h)) * sqrt(*error_face_jump_ptr);
1411  *error_face_jump_ptr = pow(*error_face_jump_ptr, 2);
1412  *error_jump_ptr += *error_face_jump_ptr;
1413  }
1414 
1415  *error_flux_ptr = 0;
1416  *error_div_ptr = 0;
1417  deltaFlux.resize(3, false);
1418  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1419  double w = getGaussPts()(3, gg) * getVolume();
1420  const double x = getCoordsAtGaussPts()(gg, 0);
1421  const double y = getCoordsAtGaussPts()(gg, 1);
1422  const double z = getCoordsAtGaussPts()(gg, 2);
1423  double flux;
1424  CHKERR cTx.getSource(fe_ent, x, y, z, flux);
1425  ;
1426  *error_div_ptr += w * pow(cTx.divergenceAtGaussPts[gg] - flux, 2);
1427  CHKERR cTx.getResistivity(fe_ent, x, y, z, invK);
1428  ;
1429  ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1430  cTx.fluxesAtGaussPts, gg);
1431  ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1433  noalias(deltaFlux) =
1434  prod(invK, flux_at_gauss_pt) + values_grad_at_gauss_pts;
1435  *error_flux_ptr += w * inner_prod(deltaFlux, deltaFlux);
1436  }
1437  *error_div_ptr = h * sqrt(*error_div_ptr);
1438  *error_div_ptr = pow(*error_div_ptr, 2);
1439  cTx.errorMap[sqrt(*error_flux_ptr + *error_div_ptr + *error_jump_ptr)] =
1440  fe_ent;
1441  // Sum/Integrate all errors
1442  cTx.sumErrorFlux += *error_flux_ptr * getVolume();
1443  cTx.sumErrorDiv += *error_div_ptr * getVolume();
1444  // FIXME: Summation should be while skeleton is calculated
1445  cTx.sumErrorJump +=
1446  *error_jump_ptr * getVolume(); /// FIXME: this need to be fixed
1448  }
1449  };
1450 
1451  /**
1452  * \brief calculate jump on entities
1453  */
1454  struct OpSkeleton
1456 
1457  /**
1458  * \brief volume element to get values from adjacent tets to face
1459  */
1461 
1462  /** store values at integration point, key of the map is sense of face in
1463  * respect to adjacent tetrahedra
1464  */
1465  map<int, VectorDouble> valMap;
1466 
1467  /**
1468  * \brief calculate values on adjacent tetrahedra to face
1469  */
1470  struct OpVolSide
1472  map<int, VectorDouble> &valMap;
1473  OpVolSide(map<int, VectorDouble> &val_map)
1475  "VALUES", UserDataOperator::OPROW),
1476  valMap(val_map) {}
1477  MoFEMErrorCode doWork(int side, EntityType type,
1480  if (data.getFieldData().size() == 0)
1482  int nb_gauss_pts = data.getN().size1();
1483  valMap[getFaceSense()].resize(nb_gauss_pts);
1484  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1485  valMap[getFaceSense()][gg] =
1486  inner_prod(trans(data.getN(gg)), data.getFieldData());
1487  }
1489  }
1490  };
1491 
1493 
1494  OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
1496  flux_name, UserDataOperator::OPROW),
1497  volSideFe(ctx.mField), cTx(ctx) {
1498  volSideFe.getOpPtrVector().push_back(new OpSkeleton::OpVolSide(valMap));
1499  }
1500 
1501  MoFEMErrorCode doWork(int side, EntityType type,
1504 
1505  if (type == MBTRI) {
1506  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1507 
1508  double def_val = 0;
1509  Tag th_error_jump;
1510  CHKERR cTx.mField.get_moab().tag_get_handle(
1511  "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1512  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1513  double *error_jump_ptr;
1514  CHKERR cTx.mField.get_moab().tag_get_by_ptr(
1515  th_error_jump, &fe_ent, 1, (const void **)&error_jump_ptr);
1516  *error_jump_ptr = 0;
1517 
1518  // check if this is essential boundary condition
1519  EntityHandle essential_bc_meshset =
1520  cTx.mField.get_finite_element_meshset("MIX_BCFLUX");
1521  if (cTx.mField.get_moab().contains_entities(essential_bc_meshset,
1522  &fe_ent, 1)) {
1523  // essential bc, np jump then, exit and go to next face
1525  }
1526 
1527  // calculate values form adjacent tets
1528  valMap.clear();
1529  CHKERR loopSideVolumes("MIX", volSideFe);
1530  ;
1531 
1532  int nb_gauss_pts = data.getN().size1();
1533 
1534  // it is only one face, so it has to be bc natural boundary condition
1535  if (valMap.size() == 1) {
1536  if (static_cast<int>(valMap.begin()->second.size()) != nb_gauss_pts) {
1537  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1538  "wrong number of integration points");
1539  }
1540  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1541  const double x = getCoordsAtGaussPts()(gg, 0);
1542  const double y = getCoordsAtGaussPts()(gg, 1);
1543  const double z = getCoordsAtGaussPts()(gg, 2);
1544  double value;
1545  CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
1546  double w = getGaussPts()(2, gg);
1547  if (static_cast<int>(getNormalsAtGaussPts().size1()) ==
1548  nb_gauss_pts) {
1549  w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1550  } else {
1551  w *= getArea();
1552  }
1553  *error_jump_ptr += w * pow(value - valMap.begin()->second[gg], 2);
1554  }
1555  } else if (valMap.size() == 2) {
1556  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1557  double w = getGaussPts()(2, gg);
1558  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1559  w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1560  } else {
1561  w *= getArea();
1562  }
1563  double delta = valMap.at(1)[gg] - valMap.at(-1)[gg];
1564  *error_jump_ptr += w * pow(delta, 2);
1565  }
1566  } else {
1567  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1568  "data inconsistency, wrong number of neighbors "
1569  "valMap.size() = %d",
1570  valMap.size());
1571  }
1572  }
1573 
1575  }
1576  };
1577 };
1578 
1579 } // namespace MixTransport
1580 
1581 #endif //_MIX_TRANPORT_ELEMENT_HPP_
1582 
1583 /***************************************************************************/ /**
1584  * \defgroup mofem_mix_transport_elem Mix transport element
1585  * \ingroup user_modules
1586  ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static const double eps
EntitiesFieldData::EntData EntData
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
@ COL
Definition: definitions.h:136
@ ROW
Definition: definitions.h:136
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:174
@ BLOCKSET
Definition: definitions.h:161
@ HVEC0
Definition: definitions.h:199
@ HVEC1
Definition: definitions.h:199
@ HVEC2
Definition: definitions.h:199
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#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
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual EntityHandle get_finite_element_meshset(const std::string &name) const =0
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreOnSideSwitch< 0 > VolumeElementForcesAndSourcesCoreOnSide
Volume element used to integrate on skeleton.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
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
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:116
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
double w(const double g, const double t)
double flux
impulse magnitude
double h
convective heat coefficient
static constexpr double delta
data for evaluation of het conductivity and heat capacity elements
Range tEts
constatins elements in block set
OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row, string val_name_col, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpError(MixTransportElement &ctx, const std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpFluxDivergenceAtGaussPts(MixTransportElement &ctx, const std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpPostProc(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
calculate values on adjacent tetrahedra to face
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
VolumeElementForcesAndSourcesCoreOnSide volSideFe
volume element to get values from adjacent tets to face
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
OpTauDotSigma_HdivHdiv(MixTransportElement &ctx, const std::string flux_name, Mat aij, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Assemble matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpVDivSigma_L2Hdiv(MixTransportElement &ctx, const std::string val_name_row, string flux_name_col, Mat aij, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
Calculate gradients of values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesGradientAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
MoFEMErrorCode postProc(const string out_file)
Post process results.
MoFEMErrorCode destroyMatrices()
destroy matrices
MoFEMErrorCode evaluateError()
Calculate error on elements.
MoFEMErrorCode solveLinearProblem()
solve problem
MoFEMErrorCode getDirichletBCIndices(IS *is)
get dof indices where essential boundary conditions are applied
MoFEMErrorCode createMatrices()
create matrices
MatrixDouble valuesGradientAtGaussPts
gradients at integration points on element
virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode calculateResidual()
calculate residual
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
VectorDouble valuesAtGaussPts
values at integration points on element
map< double, EntityHandle > errorMap
MyTriFE feTri
Instance of surface element.
virtual MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
evaluate natural (Dirichlet) boundary conditions
MyVolumeFE feVol
Instance of volume element.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
Add fields to database.
virtual MoFEMErrorCode getResistivity(const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
natural (Dirichlet) boundary conditions (set values)
MixTransportElement(MoFEM::Interface &m_field)
construction of this data structure
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
VectorDouble divergenceAtGaussPts
divergence at integration points on element
virtual MoFEMErrorCode getSource(const EntityHandle ent, const double x, const double y, const double z, double &flux)
set source term
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
MoFEMErrorCode buildProblem(BitRefLevel &ref_level)
Build problem.
virtual int get_comm_size() const =0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Mat & ts_B
Preconditioner for ts_A.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase * getVolumeFE() const
return pointer to Generic Volume Finite Element object
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Post processing.