v0.9.1
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,
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  */
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  */
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 }
Assemble constrains vector.
PetscErrorCode DMMoFEMGetSubRowIS(DM dm, IS *is)
get sub problem is
Definition: DMMMoFEM.cpp:253
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:672
Deprecated interface functions.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:288
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:706
double operator()(const double x, const double y, const double z) const
int nbIntegrationPts
number of integration points
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problemIt if true is assumed that matrix has the same indexing on rows and columns....
Definition: DMMMoFEM.cpp:378
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:221
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
int main(int argc, char *argv[])
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
const double bEta
error code
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
FTensor::Index< 'i', 3 > i
summit Index
static const bool debug
int nbCols
number if dof on column
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
MatrixDouble locMat
local entity block matrix
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
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:177
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:198
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.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
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)
Add field on domain.
Definition: Simple.cpp:324
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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)
Add field on boundary.
Definition: Simple.cpp:287
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Set integration rule to boundary elements.
bool isDiag
true if this block is on diagonal
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:212
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:457
Create finite elements instances.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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)
Add field on domain.
Definition: Simple.cpp:269
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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.
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate grad-grad operator.
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:281
constexpr int order
int nbRows
number of dofs on rows
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< ForcesAndSourcesCore > &post_proc_volume) const
Create finite element to post-process results.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:198
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:482
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
OpS(const bool beta=1)
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
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
static char help[]
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:67
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Do calculations for give operator.
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
field with C-1 continuity
Definition: definitions.h:180