v0.9.1
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
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
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
33
34  It is based on \cite arnold2006differential \cite arnold2012mixed \cite
35  reviartthomas1996
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
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
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
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  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);
282  setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
283
284  Range skeleton;
286  setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
287  moab::Interface::UNION);
289  "MIX_SKELETON");
290  }
291
292  // Define element to integrate natural boundary conditions, i.e. set values.
295  fluxes_name);
297  fluxes_name);
299  fluxes_name);
301  values_name);
302  if (mField.check_field(mesh_nodals_positions)) {
304  mesh_nodals_positions);
305  }
306
307  // Define element to apply essential boundary conditions.
310  fluxes_name);
312  fluxes_name);
314  fluxes_name);
316  values_name);
317  if (mField.check_field(mesh_nodals_positions)) {
319  mesh_nodals_positions);
320  }
321
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  */
332  // build field
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
367  // Define problem
369  // set refinement level for problem
371  // Add element to problem
374  // Boundary conditions
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");
389  }
390
391  struct OpPostProc
394  std::vector<EntityHandle> &mapGaussPts;
395  OpPostProc(moab::Interface &post_proc_mesh,
396  std::vector<EntityHandle> &map_gauss_pts)
398  "VALUES", UserDataOperator::OPCOL),
399  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
400  MoFEMErrorCode doWork(int side, EntityType type,
403  if (type != MBTET)
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  }
460  }
461  };
462
463  /**
464  * \brief Post process results
465  * @return error code
466  */
467  MoFEMErrorCode postProc(const string out_file) {
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());
480  }
481
482  Vec D, D0, F;
483  Mat Aij;
484
485  /// \brief create matrices
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);
494  }
495
496  /**
497  * \brief solve problem
498  * @return error code
499  */
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);
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
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
620  }
621
622  /// \brief calculate residual
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);
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  }
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  */
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,
690  feVol.getOpPtrVector().clear();
691  feVol.getOpPtrVector().push_back(
692  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
693  feVol.getOpPtrVector().push_back(
695  feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
696  CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
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,
706  }
707
708  /// \brief destroy matrices
711  CHKERR MatDestroy(&Aij);
712  ;
713  CHKERR VecDestroy(&D);
714  ;
715  CHKERR VecDestroy(&D0);
716  ;
717  CHKERR VecDestroy(&F);
718  ;
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  */
730
732  Mat Aij;
733  Vec F;
734
736  const std::string flux_name, Mat aij, Vec f)
738  flux_name, flux_name,
740  cTx(ctx), Aij(aij), F(f) {
741  sYmm = true;
742  }
743
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,
763  if (Aij == PETSC_NULL)
765  if (row_data.getIndices().size() == 0)
767  if (col_data.getIndices().size() == 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();
776  invK.resize(3, 3, false);
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>();
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++) {
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,
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),
821  }
822
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,
836  if (F == PETSC_NULL)
838  int nb_row = data.getIndices().size();
839  if (nb_row == 0)
841
842  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
843  Nf.resize(nb_row);
844  Nf.clear();
847  invK.resize(3, 3, false);
848  Nf.resize(nb_row);
849  Nf.clear();
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
876  }
877  };
878
879  /** \brief Assemble \f$\int_\mathcal{T} u \textrm{div}[\boldsymbol\tau] 880 * \textrm{d}\mathcal{T} \f$
881  */
884
886  Vec F;
887
888  OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row,
889  string val_name_col, Vec f)
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
900
901  MoFEMErrorCode doWork(int side, EntityType type,
904  if (data.getFieldData().size() == 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  }
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
929  }
930  };
931
932  /** \brief \f$\int_\mathcal{T} \textrm{div}[\boldsymbol\sigma] v 933 * \textrm{d}\mathcal{T} \f$
934  */
937
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)
948  val_name_row, flux_name_col,
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
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,
978  if (Aij == PETSC_NULL)
980  if (row_data.getFieldData().size() == 0)
982  if (col_data.getFieldData().size() == 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,
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),
1008  }
1009
1010  MoFEMErrorCode doWork(int side, EntityType type,
1013  if (data.getIndices().size() == 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  ;
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
1042  Vec F;
1043
1044  OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
1046  val_name, UserDataOperator::OPROW),
1047  cTx(ctx), F(f) {}
1048
1050  MoFEMErrorCode doWork(int side, EntityType type,
1053  if (data.getFieldData().size() == 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  ;
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  */
1096
1098  Vec F;
1099
1100  /**
1101  * \brief Constructor
1102  */
1103  OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name,
1104  Vec f)
1106  fluxes_name, UserDataOperator::OPROW),
1107  cTx(ctx), F(f) {}
1108
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,
1121  if (data.getFieldData().size() == 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],
1153
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  */
1174  OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
1176  flux_name, UserDataOperator::OPROW),
1177  cTx(ctx) {}
1178
1182
1183  MoFEMErrorCode doWork(int side, EntityType type,
1186  if (data.getFieldData().size() == 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
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
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
1282  }
1283  };
1284
1285  /**
1286  * \brief Calculate values at integration points
1287  */
1290
1292
1294  const std::string val_name = "VALUES")
1296  val_name, UserDataOperator::OPROW),
1297  cTx(ctx) {}
1298
1299  virtual ~OpValuesAtGaussPts() {}
1300
1301  MoFEMErrorCode doWork(int side, EntityType type,
1304  if (data.getFieldData().size() == 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
1315  }
1316  };
1317
1318  /**
1319  * \brief Calculate gradients of values at integration points
1320  */
1323
1325
1327  const std::string val_name = "VALUES")
1329  val_name, UserDataOperator::OPROW),
1330  cTx(ctx) {}
1332
1333  MoFEMErrorCode doWork(int side, EntityType type,
1336  if (data.getFieldData().size() == 0)
1338  int nb_gauss_pts = data.getDiffN().size1();
1340  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1344  prod(trans(data.getDiffN(gg)), data.getFieldData());
1345  }
1347  }
1348  };
1349
1350  /**
1351  * \brief calculate flux at integration point
1352  */
1355
1357
1359  const std::string field_name)
1361  field_name, UserDataOperator::OPROW),
1362  cTx(ctx) {}
1363
1365  MoFEMErrorCode doWork(int side, EntityType type,
1368
1369  if (data.getFieldData().size() == 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++) {
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  }
1389  }
1390  };
1391
1392  map<double, EntityHandle> errorMap;
1394  double sumErrorDiv;
1396
1397  /** \brief calculate error evaluator
1398  */
1399  struct OpError
1401
1403
1404  OpError(MixTransportElement &ctx, const std::string field_name)
1406  field_name, UserDataOperator::OPROW),
1407  cTx(ctx) {}
1408  virtual ~OpError() {}
1409
1412
1413  MoFEMErrorCode doWork(int side, EntityType type,
1416  if (type != MBTET)
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);
1488  noalias(deltaFlux) =
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
1503  }
1504  };
1505
1506  /**
1507  * \brief calculate jump on entities
1508  */
1509  struct OpSkeleton
1511
1512  /**
1513  * \brief volume element to get values from adjacent tets to face
1514  */
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
1527  map<int, VectorDouble> &valMap;
1528  OpVolSide(map<int, VectorDouble> &val_map)
1530  "VALUES", UserDataOperator::OPROW),
1531  valMap(val_map) {}
1532  MoFEMErrorCode doWork(int side, EntityType type,
1535  if (data.getFieldData().size() == 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  }
1544  }
1545  };
1546
1548
1549  OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
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,
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
1580  }
1581
1582  // calculate values form adjacent tets
1583  valMap.clear();
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
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  ******************************************************************************/
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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)
MoFEMErrorCode loopSideVolumes(const string &fe_name, VolumeElementForcesAndSourcesCoreOnSideSwitch< SWITCH > &fe_method)
MatrixDouble & getHoCoordsAtGaussPts()
coordinate at Gauss points (if hierarchical approximation of element geometry)
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
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:178
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:506
virtual int get_comm_size() const =0
VolumeElementForcesAndSourcesCoreBase * getVolumeFE() const
return pointer to Generic Volume Finite Element object
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
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:74
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
Mat & ts_B
Preconditioner for ts_A.
static constexpr double delta
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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
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:513
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")
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.
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
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
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:284
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.
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)
constexpr int order
DataForcesAndSourcesCore::EntData EntData
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
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
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
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.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:151
block name is "MAT_THERMAL"
Definition: definitions.h:229
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:97
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
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:601
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
Vec & ts_F
residual vector
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:24
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.
MatrixDouble & getHoCoordsAtGaussPts()
coordinate at Gauss points (if hierarchical approximation of element geometry)
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
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
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
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:412
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:72
OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.