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