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