v0.13.2
Loading...
Searching...
No Matches
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
13
15
16#include <PoissonOperators.hpp>
17
19
20static char help[] = "...\n\n";
21static 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 */
35struct 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/**
42 * \brief Exact gradient
43 */
44struct ExactFunctionGrad {
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;
51 return grad;
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,
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
120private:
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 /**
134 * \brief Integrate grad-grad operator
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,
195 &*locMat.data().begin(), ADD_VALUES);
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,
201 &*locMat.data().begin(), ADD_VALUES);
202 }
204 }
205};
206
207int 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>(
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
299 CHKERR simple_interface->loadFile();
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.
315 CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
316 1);
317 // Add Lagrange multiplier field on body boundary
318 CHKERR simple_interface->addBoundaryField("L", H1,
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.
324 CHKERR simple_interface->addDataField("ERROR", L2,
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
341 // Get access to PETSC-MoFEM DM manager.
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);
372 CHKERR DMMoFEMAddSubFieldRow(dm_sub_KK, "U");
373 CHKERR DMMoFEMAddSubFieldCol(dm_sub_KK, "U");
374 CHKERR DMMoFEMAddElement(dm_sub_KK,
375 simple_interface->getDomainFEName().c_str());
376 CHKERR DMMoFEMAddElement(dm_sub_KK,
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);
397 CHKERR DMMoFEMAddSubFieldRow(dm_sub_LU, "L");
398 CHKERR DMMoFEMAddSubFieldCol(dm_sub_LU, "U");
399 CHKERR DMMoFEMAddElement(dm_sub_LU,
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 }
520
521 // finish work cleaning memory, getting statistics, etc.
523
524 return 0;
525}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
Definition: adol-c_atom.cpp:46
static char help[]
static const bool debug
constexpr double a
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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:221
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:465
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:424
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMoFEM.cpp:244
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMoFEM.cpp:267
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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:554
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:503
FTensor::Index< 'i', SPACE_DIM > i
double D
FTensor::Index< 'j', 3 > j
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PetscErrorCode DMMoFEMGetSubRowIS(DM dm, IS *is)
get sub problem is
Definition: DMMoFEM.cpp:299
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
constexpr AssemblyType 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:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:327
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:320
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:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:668
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:476
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:282
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:320
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< PostProcFE > &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.