v0.13.0
analytical_poisson_field_split.cpp
Go to the documentation of this file.
1 /**
2  * \file analytical_poisson_field_split.cpp
3  * \ingroup mofem_simple_interface
4  * \example analytical_poisson_field_split.cpp
5  *
6  * For more information and detailed explain of this
7  * example see \ref poisson_tut3
8  *
9  *
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 
28 #include <PoissonOperators.hpp>
29 
30 #include <AuxPoissonFunctions.hpp>
31 
32 static char help[] = "...\n\n";
33 static const bool debug = false;
34 
35 /**
36  * \brief Function
37  *
38  * This is prescribed exact function. If this function is given by polynomial
39  * order of "p" and order of approximation is "p" or higher, solution of
40  * finite element method is exact (with machine precision).
41  *
42  * \f[
43  * u = 1+x^2+y^2+z^3
44  * \f]
45  *
46  */
47 struct ExactFunction {
48  double operator()(const double x, const double y, const double z) const {
49  return 1 + x * x + y * y + z * z * z;
50  }
51 };
52 
53 /**
54  * \brief Exact gradient
55  */
56 struct ExactFunctionGrad {
57  FTensor::Tensor1<double, 3> operator()(const double x, const double y,
58  const double z) const {
60  grad(0) = 2 * x;
61  grad(1) = 2 * y;
62  grad(2) = 3 * z * z;
63  return grad;
64  }
65 };
66 
67 /**
68  * \brief Laplacian of function.
69  *
70  * This is Laplacian of \f$u\f$, it is calculated using formula
71  * \f[
72  * \nabla^2 u(x,y,z) = \nabla \cdot \nabla u
73  * \frac{\partial^2 u}{\partial x^2}+
74  * \frac{\partial^2 u}{\partial y^2}+
75  * \frac{\partial^2 u}{\partial z^2}
76  * \f]
77  *
78  */
80  double operator()(const double x, const double y, const double z) const {
81  return 4 + 6 * z;
82  }
83 };
84 
86 
87  OpS(const bool beta = 1)
89  true),
90  bEta(beta) {}
91 
92  /**
93  * \brief Do calculations for give operator
94  * @param row_side row side number (local number) of entity on element
95  * @param col_side column side number (local number) of entity on element
96  * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
97  * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
98  * @param row_data data for row
99  * @param col_data data for column
100  * @return error code
101  */
102  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
103  EntityType col_type,
104  EntitiesFieldData::EntData &row_data,
105  EntitiesFieldData::EntData &col_data) {
107  // get number of dofs on row
108  nbRows = row_data.getIndices().size();
109  // if no dofs on row, exit that work, nothing to do here
110  if (!nbRows)
112  // get number of dofs on column
113  nbCols = col_data.getIndices().size();
114  // if no dofs on Columbia, exit nothing to do here
115  if (!nbCols)
117  // get number of integration points
118  nbIntegrationPts = getGaussPts().size2();
119  // check if entity block is on matrix diagonal
120  if (row_side == col_side && row_type == col_type) {
121  isDiag = true; // yes, it is
122  } else {
123  isDiag = false;
124  }
125  // integrate local matrix for entity block
126  CHKERR iNtegrate(row_data, col_data);
127  // assemble local matrix
128  CHKERR aSsemble(row_data, col_data);
130  }
131 
132 private:
133  const double bEta;
134 
135  ///< error code
136 
137  int nbRows; ///< number of dofs on rows
138  int nbCols; ///< number if dof on column
139  int nbIntegrationPts; ///< number of integration points
140  bool isDiag; ///< true if this block is on diagonal
141 
142  FTensor::Index<'i', 3> i; ///< summit Index
143  MatrixDouble locMat; ///< local entity block matrix
144 
145  /**
146  * \brief Integrate grad-grad operator
147  * @param row_data row data (consist base functions on row entity)
148  * @param col_data column data (consist base functions on column entity)
149  * @return error code
150  */
152  EntitiesFieldData::EntData &col_data) {
154  // set size of local entity bock
155  locMat.resize(nbRows, nbCols, false);
156  // clear matrix
157  locMat.clear();
158  // get element area
159  double area = getArea() * bEta;
160  // get integration weights
161  auto t_w = getFTensor0IntegrationWeight();
162  // get base function gradient on rows
163  auto t_row_base = row_data.getFTensor0N();
164  // loop over integration points
165  for (int gg = 0; gg != nbIntegrationPts; gg++) {
166  // take into account Jacobean
167  const double alpha = t_w * area;
168  // take fist element to local matrix
170  &*locMat.data().begin());
171  // loop over rows base functions
172  for (int rr = 0; rr != nbRows; rr++) {
173  // get column base functions gradient at gauss point gg
174  auto t_col_base = col_data.getFTensor0N(gg, 0);
175  // loop over columns
176  for (int cc = 0; cc != nbCols; cc++) {
177  // calculate element of local matrix
178  a += alpha * t_row_base * t_col_base;
179  ++t_col_base; // move to another gradient of base function on column
180  ++a; // move to another element of local matrix in column
181  }
182  ++t_row_base; // move to another element of gradient of base function on
183  // row
184  }
185  ++t_w; // move to another integration weight
186  }
188  }
189 
190  /**
191  * \brief Assemble local entity block matrix
192  * @param row_data row data (consist base functions on row entity)
193  * @param col_data column data (consist base functions on column entity)
194  * @return error code
195  */
197  EntitiesFieldData::EntData &col_data) {
199  // get pointer to first global index on row
200  const int *row_indices = &*row_data.getIndices().data().begin();
201  // get pointer to first global index on column
202  const int *col_indices = &*col_data.getIndices().data().begin();
203  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
204  : getFEMethod()->snes_B;
205  // assemble local matrix
206  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
207  &*locMat.data().begin(), ADD_VALUES);
208  if (!isDiag) {
209  // if not diagonal term and since global matrix is symmetric assemble
210  // transpose term.
211  locMat = trans(locMat);
212  CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
213  &*locMat.data().begin(), ADD_VALUES);
214  }
216  }
217 };
218 
219 int main(int argc, char *argv[]) {
220 
221  // Initialize PETSc
222  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
223  // Create MoAB database
224  moab::Core moab_core; // create database
225  moab::Interface &moab = moab_core; // create interface to database
226 
227  try {
228 
229  // Get command line options
230  int order = 3; // default approximation order
231  PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
232  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
233  "none");
234  // Set approximation order
235  CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
236  PETSC_NULL);
237  // Set testing (used by CTest)
238  CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
239  &flg_test, PETSC_NULL);
240  ierr = PetscOptionsEnd();
241  CHKERRG(ierr);
242 
243  // Create MoFEM database and link it to MoAB
244  MoFEM::Core mofem_core(moab); // create database
245  MoFEM::Interface &m_field = mofem_core; // create interface to database
246  // Register DM Manager
247  CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
248 
249  // Create vector to store approximation global error
250  Vec global_error;
251  CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
252 
253  // First we crate elements, implementation of elements is problem
254  // independent, we do not know yet what fields are present in the problem,
255  // or even we do not decided yet what approximation base or spaces we are
256  // going to use. Implementation of element is free from those constrains and
257  // can be used in different context.
258 
259  // Elements used by KSP & DM to assemble system of equations
260  boost::shared_ptr<ForcesAndSourcesCore>
261  domain_lhs_fe; ///< Volume element for the matrix
262  boost::shared_ptr<ForcesAndSourcesCore>
263  boundary_lhs_fe; ///< Boundary element for the matrix
264  boost::shared_ptr<ForcesAndSourcesCore>
265  domain_rhs_fe; ///< Volume element to assemble vector
266  boost::shared_ptr<ForcesAndSourcesCore>
267  boundary_rhs_fe; ///< Volume element to assemble vector
268  boost::shared_ptr<ForcesAndSourcesCore>
269  domain_error; ///< Volume element evaluate error
270  boost::shared_ptr<ForcesAndSourcesCore>
271  post_proc_volume; ///< Volume element to Post-process results
272  boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
273  boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
274  {
275  // Add problem specific operators the generic finite elements to calculate
276  // matrices and vectors.
279  ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
280  boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe, false);
281  // Add problem specific operators the generic finite elements to calculate
282  // error on elements and global error in H1 norm
285  global_error, domain_error);
286  // Post-process results
288  .creatFEToPostProcessResults(post_proc_volume);
289 
290  const double beta = 1;
291  boundary_penalty_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
292  new FaceElementForcesAndSourcesCore(m_field));
293  boundary_penalty_lhs_fe->getRuleHook = PoissonExample::FaceRule();
294  boundary_penalty_lhs_fe->getOpPtrVector().push_back(new OpS(beta));
295  boundary_rhs_fe->getOpPtrVector().push_back(
296  new PoissonExample::Op_g(ExactFunction(), "U", beta));
297  }
298 
299  // Get simple interface is simplified version enabling quick and
300  // easy construction of problem.
301  Simple *simple_interface;
302  // Query interface and get pointer to Simple interface
303  CHKERR m_field.getInterface(simple_interface);
304 
305  // Build problem with simple interface
306  {
307 
308  // Get options for simple interface from command line
309  CHKERR simple_interface->getOptions();
310  // Load mesh file to database
311  CHKERR simple_interface->loadFile();
312 
313  // Add field on domain and boundary. Field is declared by space and base
314  // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
315  // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
316  // number of coefficients for dof.
317  //
318  // Simple interface assumes that there are four types of field; 1) defined
319  // on body domain, 2) fields defined on body boundary, 3) skeleton field
320  // defined on finite element skeleton. Finally data field is defined on
321  // body domain. Data field is not solved but used for post-process or to
322  // keep material parameters, etc. Here we using it to calculate
323  // approximation error on elements.
324 
325  // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
326  // used to construct base functions.
327  CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
328  1);
329  // Add Lagrange multiplier field on body boundary
330  CHKERR simple_interface->addBoundaryField("L", H1,
332  // Add error (data) field, we need only L2 norm. Later order is set to 0,
333  // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
334  // element. You can use more DOFs and collate moments of error to drive
335  // anisotropic h/p-adaptivity, however this is beyond this example.
336  CHKERR simple_interface->addDataField("ERROR", L2,
338 
339  // Set fields order domain and boundary fields.
340  CHKERR simple_interface->setFieldOrder("U",
341  order); // to approximate function
342  CHKERR simple_interface->setFieldOrder("L",
343  order); // to Lagrange multipliers
344  CHKERR simple_interface->setFieldOrder(
345  "ERROR", 0); // approximation order for error
346 
347  // Setup problem. At that point database is constructed, i.e. fields,
348  // finite elements and problem data structures with local and global
349  // indexing.
350  CHKERR simple_interface->setUp();
351  }
352 
353  // Get access to PETSC-MoFEM DM manager.
354  // or more derails see
355  // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
356  // Form that point internal MoFEM data structures are linked with PETSc
357  // interface. If DM functions contains string MoFEM is is MoFEM specific DM
358  // interface function, otherwise DM function part of standard PETSc
359  // interface.
360 
361  DM dm;
362  // Get dm
363  CHKERR simple_interface->getDM(&dm);
364 
365  // Solve problem, only PETEc interface here.
366  {
367 
368  // Create the right hand side vector and vector of unknowns
369  Vec F, D;
370  CHKERR DMCreateGlobalVector(dm, &F);
371  // Create unknown vector by creating duplicate copy of F vector. only
372  // structure is duplicated no values.
373  CHKERR VecDuplicate(F, &D);
374 
375  DM dm_sub_KK, dm_sub_LU;
376  ublas::matrix<Mat> nested_matrices(2, 2);
377  ublas::vector<IS> nested_is(2);
378 
379  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_KK);
380  CHKERR DMSetType(dm_sub_KK, "DMMOFEM");
381  CHKERR DMMoFEMCreateSubDM(dm_sub_KK, dm, "SUB_KK");
382  CHKERR DMSetFromOptions(dm_sub_KK);
383  CHKERR DMMoFEMSetSquareProblem(dm_sub_KK, PETSC_TRUE);
384  CHKERR DMMoFEMAddSubFieldRow(dm_sub_KK, "U");
385  CHKERR DMMoFEMAddSubFieldCol(dm_sub_KK, "U");
386  CHKERR DMMoFEMAddElement(dm_sub_KK,
387  simple_interface->getDomainFEName().c_str());
388  CHKERR DMMoFEMAddElement(dm_sub_KK,
389  simple_interface->getBoundaryFEName().c_str());
390  CHKERR DMSetUp(dm_sub_KK);
391  CHKERR DMMoFEMGetSubRowIS(dm_sub_KK, &nested_is[0]);
392  CHKERR DMCreateMatrix(dm_sub_KK, &nested_matrices(0, 0));
393  domain_lhs_fe->ksp_B = nested_matrices(0, 0);
395  dm_sub_KK, simple_interface->getDomainFEName(), domain_lhs_fe);
396  boundary_penalty_lhs_fe->ksp_B = nested_matrices(0, 0);
398  simple_interface->getBoundaryFEName(),
399  boundary_penalty_lhs_fe);
400  CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
401  CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
402  CHKERR DMDestroy(&dm_sub_KK);
403  //
404  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_LU);
405  CHKERR DMSetType(dm_sub_LU, "DMMOFEM");
406  CHKERR DMMoFEMCreateSubDM(dm_sub_LU, dm, "SUB_LU");
407  CHKERR DMSetFromOptions(dm_sub_LU);
408  CHKERR DMMoFEMSetSquareProblem(dm_sub_LU, PETSC_FALSE);
409  CHKERR DMMoFEMAddSubFieldRow(dm_sub_LU, "L");
410  CHKERR DMMoFEMAddSubFieldCol(dm_sub_LU, "U");
411  CHKERR DMMoFEMAddElement(dm_sub_LU,
412  simple_interface->getBoundaryFEName().c_str());
413  CHKERR DMSetUp(dm_sub_LU);
414  CHKERR DMMoFEMGetSubRowIS(dm_sub_LU, &nested_is[1]);
415  CHKERR DMCreateMatrix(dm_sub_LU, &nested_matrices(1, 0));
416  boundary_lhs_fe->ksp_B = nested_matrices(1, 0);
418  dm_sub_LU, simple_interface->getBoundaryFEName(), boundary_lhs_fe);
419  CHKERR MatAssemblyBegin(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
420  CHKERR MatAssemblyEnd(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
421  // CHKERR MatCreateTranspose(nested_matrices(1,0),&nested_matrices(0,1));
422  CHKERR MatTranspose(nested_matrices(1, 0), MAT_INITIAL_MATRIX,
423  &nested_matrices(0, 1));
424  CHKERR DMDestroy(&dm_sub_LU);
425 
426  domain_rhs_fe->ksp_f = F;
427  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
428  domain_rhs_fe);
429  boundary_rhs_fe->ksp_f = F;
430  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
431  boundary_rhs_fe);
432  CHKERR VecAssemblyBegin(F);
433  CHKERR VecAssemblyEnd(F);
434 
435  Mat A;
436  nested_matrices(1, 1) = PETSC_NULL;
437 
438  if (debug) {
439  MatType type;
440  MatGetType(nested_matrices(0, 0), &type);
441  cerr << "K " << type << endl;
442  MatGetType(nested_matrices(1, 0), &type);
443  cerr << "C " << type << endl;
444  MatGetType(nested_matrices(0, 1), &type);
445  cerr << "CT " << type << endl;
446  std::string wait;
447  cerr << "UU" << endl;
448  MatView(nested_matrices(0, 0), PETSC_VIEWER_DRAW_WORLD);
449  std::cin >> wait;
450  cerr << "LU" << endl;
451  MatView(nested_matrices(1, 0), PETSC_VIEWER_DRAW_WORLD);
452  std::cin >> wait;
453  cerr << "UL" << endl;
454  MatView(nested_matrices(0, 1), PETSC_VIEWER_DRAW_WORLD);
455  std::cin >> wait;
456  }
457 
458  CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is[0], 2, &nested_is[0],
459  &nested_matrices(0, 0), &A);
460 
461  // Create solver and link it to DM
462  KSP solver;
463  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
464  CHKERR KSPSetFromOptions(solver);
465  // Set operators
466  CHKERR KSPSetOperators(solver, A, A);
467  PC pc;
468  CHKERR KSPGetPC(solver, &pc);
469  PetscBool is_pcfs = PETSC_FALSE;
470  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
471  if (is_pcfs) {
472  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[0]);
473  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[1]);
474  } else {
475  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
476  "Only works with pre-conditioner PCFIELDSPLIT");
477  }
478  // Set-up solver, is type of solver and pre-conditioners
479  CHKERR KSPSetUp(solver);
480  // At solution process, KSP solver using DM creates matrices, Calculate
481  // values of the left hand side and the right hand side vector. then
482  // solves system of equations. Results are stored in vector D.
483  CHKERR KSPSolve(solver, F, D);
484 
485  // Scatter solution on the mesh. Stores unknown vector on field on the
486  // mesh.
487  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
488 
489  // Clean data. Solver and vector are not needed any more.
490  CHKERR KSPDestroy(&solver);
491  for (int i = 0; i != 2; i++) {
492  CHKERR ISDestroy(&nested_is[i]);
493  for (int j = 0; j != 2; j++) {
494  CHKERR MatDestroy(&nested_matrices(i, j));
495  }
496  }
497  CHKERR MatDestroy(&A);
498  CHKERR VecDestroy(&D);
499  CHKERR VecDestroy(&F);
500  }
501 
502  // Calculate error
503  {
504  // Loop over all elements in mesh, and run users operators on each
505  // element.
506  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
507  domain_error);
509  global_error);
510  CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
511  if (flg_test == PETSC_TRUE) {
512  CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
513  }
514  }
515 
516  {
517  // Loop over all elements in the mesh and for each execute
518  // post_proc_volume element and operators on it.
519  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
520  post_proc_volume);
521  // Write results
522  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
523  post_proc_volume)
524  ->writeFile("out_vol.h5m");
525  }
526 
527  // Destroy DM, no longer needed.
528  CHKERR DMDestroy(&dm);
529 
530  // Destroy ghost vector
531  CHKERR VecDestroy(&global_error);
532  }
533  CATCH_ERRORS;
534 
535  // finish work cleaning memory, getting statistics, etc.
537 
538  return 0;
539 }
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main(int argc, char *argv[])
static char help[]
static const bool debug
constexpr double a
EntitiesFieldData::EntData EntData
constexpr double alpha
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#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
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMMoFEM.cpp:213
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMMoFEM.cpp:414
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:234
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:257
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:493
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
PetscErrorCode DMMoFEMGetSubRowIS(DM dm, IS *is)
get sub problem is
Definition: DMMMoFEM.cpp:289
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
double A
Exact gradient.
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double z) const
double operator()(const double x, const double y, const double z) const
double operator()(const double x, const double y, const double z) const
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
int nbIntegrationPts
number of integration points
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
OpS(const bool beta=1)
int nbRows
number of dofs on rows
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
int nbCols
number if dof on column
const double bEta
error code
MatrixDouble locMat
local entity block matrix
bool isDiag
true if this block is on diagonal
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
FTensor::Index< 'i', 3 > i
summit Index
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
Create finite elements instances.
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< ForcesAndSourcesCore > &post_proc_volume) const
Create finite element to post-process results.
MoFEMErrorCode createFEToEvaluateError(boost::function< double(const double, const double, const double)> f_u, boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> g_u, Vec global_error, boost::shared_ptr< ForcesAndSourcesCore > &domain_error) const
Create finite element to calculate error.
MoFEMErrorCode createFEToAssembleMatrixAndVector(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, bool trans=true) const
Create finite element to calculate matrix and vectors.
Set integration rule to boundary elements.
Assemble constrains vector.