v0.13.0
h_adaptive_transport.cpp
Go to the documentation of this file.
1 /** \file h_adaptive_transport.cpp
2 \brief Example implementation of transport problem using mixed formulation
3 
4 \todo Should be implemented and tested problem from this article
5 Demkowicz, Leszek, and Jayadeep Gopalakrishnan. "Analysis of the DPG method for
6 the Poisson equation." SIAM Journal on Numerical Analysis 49.5 (2011):
7 1788-1809.
8 
9 \ingroup mofem_mix_transport_elem
10 */
11 
12 /* This file is part of MoFEM.
13  * MoFEM is free software: you can redistribute it and/or modify it under
14  * the terms of the GNU Lesser General Public License as published by the
15  * Free Software Foundation, either version 3 of the License, or (at your
16  * option) any later version.
17  *
18  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
19  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21  * License for more details.
22  *
23  * You should have received a copy of the GNU Lesser General Public
24  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
25 
26 #include <BasicFiniteElements.hpp>
27 #include <MixTransportElement.hpp>
28 
29 using namespace MoFEM;
30 using namespace MixTransport;
31 
32 static char help[] = "-my_file input file"
33  "-my_order order of approximation"
34  "-nb_levels number of refinements levels"
35  "\n\n";
36 
37 /**
38  * Data structure to pass information between function evaluating boundary
39  * values and fluxes and generic data structures for boundary conditions on
40  * meshsets.
41  */
42 struct BcFluxData {
43  Range eNts;
44  double fLux;
45 };
46 typedef map<int, BcFluxData> BcFluxMap;
47 
48 /** \brief Application of mix transport data structure
49  *
50  * MixTransportElement is a class collecting functions, operators and
51  * data for mix implementation of transport element. See there to
52  * learn how elements are created or how operators look like.
53  *
54  * Some methods in MixTransportElement are abstract, f.e. user need to
55  * implement own source therm.
56 
57  * \ingroup mofem_mix_transport_elem
58  */
59 struct MyTransport : public MixTransportElement {
60 
63  double lastFlux;
64 
65  MyTransport(MoFEM::Interface &m_field, BcFluxMap &bc_flux_map)
66  : MixTransportElement(m_field), bcFluxMap(bc_flux_map), lastEnt(0),
67  lastFlux(0) {}
68 
69  /**
70  * \brief set source term
71  * @param ent handle to entity on which function is evaluated
72  * @param x coord
73  * @param y coord
74  * @param z coord
75  * @param flux reference to source term set by function
76  * @return error code
77  */
78  MoFEMErrorCode getSource(EntityHandle ent, const double x, const double y,
79  const double z, double &flux) {
81  flux = 0;
83  }
84 
85  /**
86  * \brief natural (Dirihlet) boundary conditions (set values)
87  * @param ent handle to finite element entity
88  * @param x coord
89  * @param y coord
90  * @param z coord
91  * @param value reference to value set by function
92  * @return error code
93  */
94  MoFEMErrorCode getBcOnValues(const EntityHandle ent, const double x,
95  const double y, const double z, double &value) {
97  value = 0;
99  }
100 
101  /**
102  * \brief essential (Neumann) boundary condition (set fluxes)
103  * @param ent handle to finite element entity
104  * @param x coord
105  * @param y coord
106  * @param z coord
107  * @param flux reference to flux which is set by function
108  * @return [description]
109  */
110  MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
111  const double y, const double z, double &flux) {
113  if (lastEnt == ent) {
114  flux = lastFlux;
115  } else {
116  flux = 0;
117  for (BcFluxMap::iterator mit = bcFluxMap.begin(); mit != bcFluxMap.end();
118  mit++) {
119  Range &tris = mit->second.eNts;
120  if (tris.find(ent) != tris.end()) {
121  flux = mit->second.fLux;
122  }
123  }
124  lastEnt = ent;
125  lastFlux = flux;
126  }
128  }
129 
130  /**
131  * \brief set-up boundary conditions
132  * @param ref_level mesh refinement level
133 
134  \note It is assumed that user would like to something non-standard with
135  boundary conditions, have a own type of data structures to pass to functions
136  calculating values and fluxes on boundary. For example BcFluxMap. That way
137  this function is implemented here not in generic class MixTransportElement.
138 
139  * @return error code
140  */
143  Range tets;
144  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
145  ref_level, BitRefLevel().set(), MBTET, tets);
146  Skinner skin(&mField.get_moab());
147  Range skin_faces; // skin faces from 3d ents
148  CHKERR skin.find_skin(0, tets, false, skin_faces);
149  // note: what is essential (dirichlet) is natural (neumann) for mix-FE
150  // compared to classical FE
151  Range natural_bc;
153  mField, NODESET | TEMPERATURESET, it)) {
154  Range tris;
155  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, tris,
156  true);
157  natural_bc.insert(tris.begin(), tris.end());
158  }
160  mField, SIDESET | HEATFLUXSET, it)) {
161  HeatFluxCubitBcData mydata;
162  CHKERR it->getBcDataStructure(mydata);
163  if (mydata.data.flag1 == 1) {
164  Range tris;
165  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, tris,
166  true);
167  bcFluxMap[it->getMeshsetId()].eNts = tris;
168  bcFluxMap[it->getMeshsetId()].fLux = mydata.data.value1;
169  // cerr << bcFluxMap[it->getMeshsetId()].eNts << endl;
170  // cerr << bcFluxMap[it->getMeshsetId()].fLux << endl;
171  }
172  }
173  Range essential_bc = subtract(skin_faces, natural_bc);
174  Range bit_tris;
175  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
176  ref_level, BitRefLevel().set(), MBTRI, bit_tris);
177  essential_bc = intersect(bit_tris, essential_bc);
178  natural_bc = intersect(bit_tris, natural_bc);
179  CHKERR mField.add_ents_to_finite_element_by_type(essential_bc, MBTRI,
180  "MIX_BCFLUX");
181  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc, MBTRI,
182  "MIX_BCVALUE");
183  // CHKERR
184  // mField.add_ents_to_finite_element_by_type(skin_faces,MBTRI,"MIX_BCVALUE");
186  }
187 
188  /**
189  * \brief Refine mesh
190  * @param ufe general data structure
191  * @param nb_levels number of refinement levels
192  * @param order set order of approximation
193  * @return errpr code
194 
195  Refinement of could result in distorted mesh, for example, imagine when you
196  have two levels of non-uniform refinement. Some tetrahedra on the mesh at
197  first refinement instance are only refined by splitting subset of edges on
198  it. Then refined child tetrahedra usually will have worse quality than
199  quality of parent element. Refining such element in subsequent mesh
200  refinement, potentially will deteriorate elements quality even worse. To
201  prevent that adding new refinement level, recreate whole hierarchy of meshes.
202 
203  Note on subsequent improvement could include refinement of
204  tetrahedra from different levels, including initial mesh. So refinement two
205  could split elements created during refinement one and also split elements
206  from an initial mesh.
207 
208  That adding the new refinement level creates refinement hierarchy of meshes
209  from a scratch, not adding to existing one.
210 
211  Entities from previous hierarchy are used in that process, but bit levels on
212  those entities are squashed.
213 
214  */
215  MoFEMErrorCode refineMesh(MixTransportElement &ufe, const int nb_levels,
216  const int order) {
217  MeshRefinement *refine_ptr;
219  // get refined edges having child vertex
220  auto ref_ents_ptr = mField.get_ref_ents();
223  const RefEntsByComposite &ref_ents =
224  ref_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
225  RefEntsByComposite::iterator rit, hi_rit;
226  rit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
227  hi_rit = ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
228  Range refined_edges;
229  // thist loop is over vertices which parent is edge
230  for (; rit != hi_rit; rit++) {
231  refined_edges.insert((*rit)->getParentEnt()); // get parent edge
232  }
233  // get tets which has large error
234  Range tets_to_refine;
235  const double max_error = ufe.errorMap.rbegin()->first;
236  // int size = ((double)5/6)*ufe.errorMap.size();
237  for (map<double, EntityHandle>::iterator mit = ufe.errorMap.begin();
238  mit != ufe.errorMap.end(); mit++) {
239  // cerr << mit->first << " " << mit->second << endl;
240  // if((size--)>0) continue;
241  if (mit->first < 0.25 * max_error)
242  continue;
243  tets_to_refine.insert(mit->second);
244  }
245  Range tets_to_refine_edges;
246  CHKERR mField.get_moab().get_adjacencies(
247  tets_to_refine, 1, false, tets_to_refine_edges, moab::Interface::UNION);
248  refined_edges.merge(tets_to_refine_edges);
249  CHKERR mField.getInterface(refine_ptr);
250  for (int ll = 0; ll != nb_levels; ll++) {
251  Range edges;
252  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
253  BitRefLevel().set(ll), BitRefLevel().set(), MBEDGE, edges);
254  edges = intersect(edges, refined_edges);
255  // add edges to refine at current level edges (some of the where refined
256  // before)
258  edges, BitRefLevel().set(ll + 1));
259  // get tets at current level
260  Range tets;
261  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
262  BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
263  CHKERR refine_ptr->refineTets(tets, BitRefLevel().set(ll + 1));
264  CHKERR updateMeshsetsFieldsAndElements(ll + 1);
265  }
266 
267  // update fields and elements
268  EntityHandle ref_meshset;
269  CHKERR mField.get_moab().create_meshset(MESHSET_SET, ref_meshset);
270  {
271  // cerr << BitRefLevel().set(nb_levels) << endl;
272  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
273  BitRefLevel().set(nb_levels), BitRefLevel().set(), MBTET,
274  ref_meshset);
275 
276  Range ref_tets;
277  CHKERR mField.get_moab().get_entities_by_type(ref_meshset, MBTET,
278  ref_tets);
279 
280  // add entities to field
281  CHKERR mField.add_ents_to_field_by_type(ref_meshset, MBTET, "FLUXES");
282  CHKERR mField.add_ents_to_field_by_type(ref_meshset, MBTET, "VALUES");
283  CHKERR mField.set_field_order(0, MBTET, "FLUXES", order + 1);
284  CHKERR mField.set_field_order(0, MBTRI, "FLUXES", order + 1);
285  CHKERR mField.set_field_order(0, MBTET, "VALUES", order);
286 
287  // add entities to skeleton
288  Range ref_tris;
289  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
290  BitRefLevel().set(nb_levels), BitRefLevel().set(), MBTRI, ref_tris);
291  CHKERR mField.add_ents_to_finite_element_by_type(ref_tris, MBTRI,
292  "MIX_SKELETON");
293 
294  // add entities to finite elements
296  mField, BLOCKSET | MAT_THERMALSET, it)) {
297  Mat_Thermal temp_data;
298  CHKERR it->getAttributeDataStructure(temp_data);
299  setOfBlocks[it->getMeshsetId()].cOnductivity =
300  temp_data.data.Conductivity;
301  setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
302  CHKERR mField.get_moab().get_entities_by_type(
303  it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
304  setOfBlocks[it->getMeshsetId()].tEts =
305  intersect(ref_tets, setOfBlocks[it->getMeshsetId()].tEts);
306  CHKERR mField.add_ents_to_finite_element_by_type(
307  setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
308  }
309  }
310  CHKERR mField.get_moab().delete_entities(&ref_meshset, 1);
312  }
313 
314  /**
315  * \brief Squash bits of entities
316 
317  Information about hierarchy of meshsets is lost, but entities are not deleted
318  from the mesh. After squash entities bits, new hierarchy can be created.
319 
320  * @return error code
321  */
324  BitRefLevel all_but_0;
325  all_but_0.set(0);
326  all_but_0.flip();
327  BitRefLevel garbage_bit;
328  garbage_bit.set(BITREFLEVEL_SIZE - 1); // Garbage level
329  auto ref_ents_ptr = mField.get_ref_ents();
330  RefEntity_multiIndex::iterator mit = ref_ents_ptr->begin();
331  for (; mit != ref_ents_ptr->end(); mit++) {
332  if (mit->get()->getEntType() == MBENTITYSET)
333  continue;
334  BitRefLevel bit = mit->get()->getBitRefLevel();
335  if ((all_but_0 & bit) == bit) {
336  *(const_cast<RefEntity *>(mit->get())->getBitRefLevelPtr()) =
337  garbage_bit;
338  } else {
339  *(const_cast<RefEntity *>(mit->get())->getBitRefLevelPtr()) =
340  BitRefLevel().set(0);
341  }
342  }
344  }
345 
346  /**
347  * \brief update meshsets with new entities after mesh refinement
348  * @param nb_levels nb_levels
349  * @param order appropriate order
350  * @return error code
351  */
353  BitRefLevel ref_level;
355  ref_level.set(nb_levels);
356  for (_IT_CUBITMESHSETS_FOR_LOOP_(mField, it)) {
357  EntityHandle meshset = it->meshset;
358  CHKERR mField.getInterface<BitRefManager>()
359  ->updateMeshsetByEntitiesChildren(meshset, ref_level, meshset,
360  MBVERTEX, true);
361  CHKERR mField.getInterface<BitRefManager>()
362  ->updateMeshsetByEntitiesChildren(meshset, ref_level, meshset, MBEDGE,
363  true);
364  CHKERR mField.getInterface<BitRefManager>()
365  ->updateMeshsetByEntitiesChildren(meshset, ref_level, meshset, MBTRI,
366  true);
367  CHKERR mField.getInterface<BitRefManager>()
368  ->updateMeshsetByEntitiesChildren(meshset, ref_level, meshset, MBTET,
369  true);
370  }
372  }
373 };
374 
375 int main(int argc, char *argv[]) {
376 
377  const string default_options = "-ksp_type fgmres \n"
378  "-pc_type lu \n"
379  "-pc_factor_mat_solver_type mumps \n"
380  "-mat_mumps_icntl_20 0 \n"
381  "-ksp_monitor\n";
382 
383  string param_file = "param_file.petsc";
384  if (!static_cast<bool>(ifstream(param_file))) {
385  std::ofstream file(param_file.c_str(), std::ios::ate);
386  if (file.is_open()) {
387  file << default_options;
388  file.close();
389  }
390  }
391 
392  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
393 
394  try {
395 
396  moab::Core mb_instance;
397  moab::Interface &moab = mb_instance;
398  int rank;
399  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
400 
401  // get file name form command line
402  PetscBool flg = PETSC_TRUE;
403  char mesh_file_name[255];
404  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
405  mesh_file_name, 255, &flg);
406  if (flg != PETSC_TRUE) {
407  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
408  "*** ERROR -my_file (MESH FILE NEEDED)");
409  }
410 
411  const char *option;
412  option = "";
413  CHKERR moab.load_file(mesh_file_name, 0, option);
414 
415  // Create mofem interface
416  MoFEM::Core core(moab);
417  MoFEM::Interface &m_field = core;
418 
419  // Add meshsets with material and boundary conditions
420  MeshsetsManager *meshsets_manager_ptr;
421  CHKERR m_field.getInterface(meshsets_manager_ptr);
422  CHKERR meshsets_manager_ptr->setMeshsetFromFile();
423 
424  PetscPrintf(PETSC_COMM_WORLD,
425  "Read meshsets add added meshsets for bc.cfg\n");
426  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
427  PetscPrintf(PETSC_COMM_WORLD, "%s",
428  static_cast<std::ostringstream &>(
429  std::ostringstream().seekp(0) << *it << endl)
430  .str()
431  .c_str());
432  cerr << *it << endl;
433  }
434 
435  // set entities bit level
436  BitRefLevel ref_level;
437  ref_level.set(0);
438  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
439  0, 3, ref_level);
440 
441  // set app. order
442  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
443  // (Mark Ainsworth & Joe Coyle)
444  PetscInt order;
445  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
446  &flg);
447  if (flg != PETSC_TRUE) {
448  order = 0;
449  }
450 
451  // finite elements
452 
453  BcFluxMap bc_flux_map;
454  MyTransport ufe(m_field, bc_flux_map);
455 
456  // Initially calculate problem on coarse mesh
457 
458  CHKERR ufe.addFields("VALUES", "FLUXES", order);
459  CHKERR ufe.addFiniteElements("FLUXES", "VALUES");
460  // Set boundary conditions
461  CHKERR ufe.addBoundaryElements(ref_level);
462  CHKERR ufe.buildProblem(ref_level);
463  CHKERR ufe.createMatrices();
466  CHKERR ufe.evaluateError();
467  CHKERR ufe.destroyMatrices();
468  CHKERR ufe.postProc("out_0.h5m");
469 
470  int nb_levels = 5; // default number of refinement levels
471  // get number of refinement levels form command line
472  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-nb_levels", &nb_levels,
473  PETSC_NULL);
474 
475  // refine mesh, solve problem and do it again until number of refinement
476  // levels are exceeded.
477  for (int ll = 1; ll != nb_levels; ll++) {
478  const int nb_levels = ll;
479  CHKERR ufe.squashBits();
480  CHKERR ufe.refineMesh(ufe, nb_levels, order);
481  ref_level = BitRefLevel().set(nb_levels);
482  bc_flux_map.clear();
483  CHKERR ufe.addBoundaryElements(ref_level);
484  CHKERR ufe.buildProblem(ref_level);
485  CHKERR ufe.createMatrices();
488  CHKERR ufe.evaluateError();
489  CHKERR ufe.destroyMatrices();
490  CHKERR ufe.postProc(
491  static_cast<std::ostringstream &>(std::ostringstream().seekp(0)
492  << "out_" << nb_levels << ".h5m")
493  .str());
494  }
495  }
496  CATCH_ERRORS;
497 
499 
500  return 0;
501 }
const std::string default_options
std::string param_file
Mix implementation of transport element.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ TEMPERATURESET
Definition: definitions.h:168
@ HEATFLUXSET
Definition: definitions.h:169
@ NODESET
Definition: definitions.h:159
@ SIDESET
Definition: definitions.h:160
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:174
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
int main(int argc, char *argv[])
static char help[]
map< int, BcFluxData > BcFluxMap
auto bit
set bit
char mesh_file_name[255]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
double flux
impulse magnitude
MoFEMErrorCode postProc(const string out_file)
Post process results.
MoFEMErrorCode destroyMatrices()
destroy matrices
MoFEMErrorCode evaluateError()
Calculate error on elements.
MoFEMErrorCode solveLinearProblem()
solve problem
MoFEMErrorCode createMatrices()
create matrices
MoFEMErrorCode calculateResidual()
calculate residual
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
map< double, EntityHandle > errorMap
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
Add fields to database.
MoFEMErrorCode buildProblem(BitRefLevel &ref_level)
Build problem.
Managing BitRefLevels.
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Definition of the heat flux bc data structure.
Definition: BCData.hpp:437
Thermal material data structure.
Mesh refinement interface.
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
Struct keeps handle to refined handle.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Application of mix transport data structure.
MoFEMErrorCode updateMeshsetsFieldsAndElements(const int nb_levels)
update meshsets with new entities after mesh refinement
MyTransport(MoFEM::Interface &m_field, BcFluxMap &bc_flux_map)
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const double x, const double y, const double z, double &value)
natural (Dirihlet) boundary conditions (set values)
MoFEMErrorCode addBoundaryElements(BitRefLevel &ref_level)
set-up boundary conditions
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode squashBits()
Squash bits of entities.
EntityHandle lastEnt
MoFEMErrorCode refineMesh(MixTransportElement &ufe, const int nb_levels, const int order)
Refine mesh.
MoFEMErrorCode getSource(EntityHandle ent, const double x, const double y, const double z, double &flux)
set source term
BcFluxMap & bcFluxMap