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