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