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