v0.15.0
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
14#include <MoFEM.hpp>
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
73struct OpS : public FaceElementForcesAndSourcesCore::UserDataOperator {
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_NULLPTR ? 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 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
221 "none");
222 // Set approximation order
223 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
224 PETSC_NULLPTR);
225 // Set testing (used by CTest)
226 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
227 &flg_test, PETSC_NULLPTR);
228 PetscOptionsEnd();
229
230 // Create MoFEM database and link it to MoAB
231 MoFEM::Core mofem_core(moab); // create database
232 MoFEM::Interface &m_field = mofem_core; // create interface to database
233 // Register DM Manager
234 CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
235
236 // Create vector to store approximation global error
237 Vec global_error;
238 CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
239
240 // First we crate elements, implementation of elements is problem
241 // independent, we do not know yet what fields are present in the problem,
242 // or even we do not decided yet what approximation base or spaces we are
243 // going to use. Implementation of element is free from those constrains and
244 // can be used in different context.
245
246 // Elements used by KSP & DM to assemble system of equations
247 boost::shared_ptr<ForcesAndSourcesCore>
248 domain_lhs_fe; ///< Volume element for the matrix
249 boost::shared_ptr<ForcesAndSourcesCore>
250 boundary_lhs_fe; ///< Boundary element for the matrix
251 boost::shared_ptr<ForcesAndSourcesCore>
252 domain_rhs_fe; ///< Volume element to assemble vector
253 boost::shared_ptr<ForcesAndSourcesCore>
254 boundary_rhs_fe; ///< Volume element to assemble vector
255 boost::shared_ptr<ForcesAndSourcesCore>
256 domain_error; ///< Volume element evaluate error
257 boost::shared_ptr<PoissonExample::PostProcFE>
258 post_proc_volume; ///< Volume element to Post-process results
259 boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
260 boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
261 {
262 // Add problem specific operators the generic finite elements to calculate
263 // matrices and vectors.
266 ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
267 boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe, false);
268 // Add problem specific operators the generic finite elements to calculate
269 // error on elements and global error in H1 norm
272 global_error, domain_error);
273 // Post-process results
275 .creatFEToPostProcessResults(post_proc_volume);
276
277 const double beta = 1;
278 boundary_penalty_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
280 boundary_penalty_lhs_fe->getRuleHook = PoissonExample::FaceRule();
281 boundary_penalty_lhs_fe->getOpPtrVector().push_back(new OpS(beta));
282 boundary_rhs_fe->getOpPtrVector().push_back(
283 new PoissonExample::Op_g(ExactFunction(), "U", beta));
284 }
285
286 // Get simple interface is simplified version enabling quick and
287 // easy construction of problem.
288 Simple *simple_interface;
289 // Query interface and get pointer to Simple interface
290 CHKERR m_field.getInterface(simple_interface);
291
292 // Build problem with simple interface
293 {
294
295 // Get options for simple interface from command line
296 CHKERR simple_interface->getOptions();
297 // Load mesh file to database
298 CHKERR simple_interface->loadFile();
299
300 // Add field on domain and boundary. Field is declared by space and base
301 // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
302 // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
303 // number of coefficients for dof.
304 //
305 // Simple interface assumes that there are four types of field; 1) defined
306 // on body domain, 2) fields defined on body boundary, 3) skeleton field
307 // defined on finite element skeleton. Finally data field is defined on
308 // body domain. Data field is not solved but used for post-process or to
309 // keep material parameters, etc. Here we using it to calculate
310 // approximation error on elements.
311
312 // Add domain field "U" in space H1 and Legendre base, Ainsworth recipe is
313 // used to construct base functions.
314 CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
315 1);
316 // Add Lagrange multiplier field on body boundary
317 CHKERR simple_interface->addBoundaryField("L", H1,
319 // Add error (data) field, we need only L2 norm. Later order is set to 0,
320 // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
321 // element. You can use more DOFs and collate moments of error to drive
322 // anisotropic h/p-adaptivity, however this is beyond this example.
323 CHKERR simple_interface->addDataField("ERROR", L2,
325
326 // Set fields order domain and boundary fields.
327 CHKERR simple_interface->setFieldOrder("U",
328 order); // to approximate function
329 CHKERR simple_interface->setFieldOrder("L",
330 order); // to Lagrange multipliers
331 CHKERR simple_interface->setFieldOrder(
332 "ERROR", 0); // approximation order for error
333
334 // Setup problem. At that point database is constructed, i.e. fields,
335 // finite elements and problem data structures with local and global
336 // indexing.
337 CHKERR simple_interface->setUp();
338 }
339
340 // Get access to PETSC-MoFEM DM manager.
341 // or more derails see
342 // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
343 // Form that point internal MoFEM data structures are linked with PETSc
344 // interface. If DM functions contains string MoFEM is is MoFEM specific DM
345 // interface function, otherwise DM function part of standard PETSc
346 // interface.
347
348 DM dm;
349 // Get dm
350 CHKERR simple_interface->getDM(&dm);
351
352 // Solve problem, only PETEc interface here.
353 {
354
355 // Create the right hand side vector and vector of unknowns
356 Vec F, D;
357 CHKERR DMCreateGlobalVector(dm, &F);
358 // Create unknown vector by creating duplicate copy of F vector. only
359 // structure is duplicated no values.
360 CHKERR VecDuplicate(F, &D);
361
362 DM dm_sub_KK, dm_sub_LU;
363 ublas::matrix<Mat> nested_matrices(2, 2);
364 ublas::vector<IS> nested_is(2);
365
366 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_KK);
367 CHKERR DMSetType(dm_sub_KK, "DMMOFEM");
368 CHKERR DMMoFEMCreateSubDM(dm_sub_KK, dm, "SUB_KK");
369 CHKERR DMSetFromOptions(dm_sub_KK);
370 CHKERR DMMoFEMSetSquareProblem(dm_sub_KK, PETSC_TRUE);
371 CHKERR DMMoFEMAddSubFieldRow(dm_sub_KK, "U");
372 CHKERR DMMoFEMAddSubFieldCol(dm_sub_KK, "U");
373 CHKERR DMMoFEMAddElement(dm_sub_KK,
374 simple_interface->getDomainFEName().c_str());
375 CHKERR DMMoFEMAddElement(dm_sub_KK,
376 simple_interface->getBoundaryFEName().c_str());
377 CHKERR DMSetUp(dm_sub_KK);
378 CHKERR DMMoFEMGetSubRowIS(dm_sub_KK, &nested_is[0]);
379 CHKERR DMCreateMatrix(dm_sub_KK, &nested_matrices(0, 0));
380 domain_lhs_fe->ksp_B = nested_matrices(0, 0);
382 dm_sub_KK, simple_interface->getDomainFEName(), domain_lhs_fe);
383 boundary_penalty_lhs_fe->ksp_B = nested_matrices(0, 0);
385 simple_interface->getBoundaryFEName(),
386 boundary_penalty_lhs_fe);
387 CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
388 CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
389 CHKERR DMDestroy(&dm_sub_KK);
390 //
391 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_LU);
392 CHKERR DMSetType(dm_sub_LU, "DMMOFEM");
393 CHKERR DMMoFEMCreateSubDM(dm_sub_LU, dm, "SUB_LU");
394 CHKERR DMSetFromOptions(dm_sub_LU);
395 CHKERR DMMoFEMSetSquareProblem(dm_sub_LU, PETSC_FALSE);
396 CHKERR DMMoFEMAddSubFieldRow(dm_sub_LU, "L");
397 CHKERR DMMoFEMAddSubFieldCol(dm_sub_LU, "U");
398 CHKERR DMMoFEMAddElement(dm_sub_LU,
399 simple_interface->getBoundaryFEName().c_str());
400 CHKERR DMSetUp(dm_sub_LU);
401 CHKERR DMMoFEMGetSubRowIS(dm_sub_LU, &nested_is[1]);
402 CHKERR DMCreateMatrix(dm_sub_LU, &nested_matrices(1, 0));
403 boundary_lhs_fe->ksp_B = nested_matrices(1, 0);
405 dm_sub_LU, simple_interface->getBoundaryFEName(), boundary_lhs_fe);
406 CHKERR MatAssemblyBegin(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
407 CHKERR MatAssemblyEnd(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
408 // CHKERR MatCreateTranspose(nested_matrices(1,0),&nested_matrices(0,1));
409 CHKERR MatTranspose(nested_matrices(1, 0), MAT_INITIAL_MATRIX,
410 &nested_matrices(0, 1));
411 CHKERR DMDestroy(&dm_sub_LU);
412
413 domain_rhs_fe->ksp_f = F;
414 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
415 domain_rhs_fe);
416 boundary_rhs_fe->ksp_f = F;
417 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
418 boundary_rhs_fe);
419 CHKERR VecAssemblyBegin(F);
420 CHKERR VecAssemblyEnd(F);
421
422 Mat A;
423 nested_matrices(1, 1) = PETSC_NULLPTR;
424
425 if (debug) {
426 MatType type;
427 MatGetType(nested_matrices(0, 0), &type);
428 cerr << "K " << type << endl;
429 MatGetType(nested_matrices(1, 0), &type);
430 cerr << "C " << type << endl;
431 MatGetType(nested_matrices(0, 1), &type);
432 cerr << "CT " << type << endl;
433 std::string wait;
434 cerr << "UU" << endl;
435 MatView(nested_matrices(0, 0), PETSC_VIEWER_DRAW_WORLD);
436 std::cin >> wait;
437 cerr << "LU" << endl;
438 MatView(nested_matrices(1, 0), PETSC_VIEWER_DRAW_WORLD);
439 std::cin >> wait;
440 cerr << "UL" << endl;
441 MatView(nested_matrices(0, 1), PETSC_VIEWER_DRAW_WORLD);
442 std::cin >> wait;
443 }
444
445 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is[0], 2, &nested_is[0],
446 &nested_matrices(0, 0), &A);
447
448 // Create solver and link it to DM
449 KSP solver;
450 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
451 CHKERR KSPSetFromOptions(solver);
452 // Set operators
453 CHKERR KSPSetOperators(solver, A, A);
454 PC pc;
455 CHKERR KSPGetPC(solver, &pc);
456 PetscBool is_pcfs = PETSC_FALSE;
457 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
458 if (is_pcfs) {
459 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[0]);
460 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[1]);
461 } else {
462 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
463 "Only works with pre-conditioner PCFIELDSPLIT");
464 }
465 // Set-up solver, is type of solver and pre-conditioners
466 CHKERR KSPSetUp(solver);
467 // At solution process, KSP solver using DM creates matrices, Calculate
468 // values of the left hand side and the right hand side vector. then
469 // solves system of equations. Results are stored in vector D.
470 CHKERR KSPSolve(solver, F, D);
471
472 // Scatter solution on the mesh. Stores unknown vector on field on the
473 // mesh.
474 CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
475
476 // Clean data. Solver and vector are not needed any more.
477 CHKERR KSPDestroy(&solver);
478 for (int i = 0; i != 2; i++) {
479 CHKERR ISDestroy(&nested_is[i]);
480 for (int j = 0; j != 2; j++) {
481 CHKERR MatDestroy(&nested_matrices(i, j));
482 }
483 }
484 CHKERR MatDestroy(&A);
485 CHKERR VecDestroy(&D);
486 CHKERR VecDestroy(&F);
487 }
488
489 // Calculate error
490 {
491 // Loop over all elements in mesh, and run users operators on each
492 // element.
493 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
494 domain_error);
496 global_error);
497 CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
498 if (flg_test == PETSC_TRUE) {
499 CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
500 }
501 }
502
503 {
504 // Loop over all elements in the mesh and for each execute
505 // post_proc_volume element and operators on it.
506 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
507 post_proc_volume);
508 // Write results
509 post_proc_volume->writeFile("out_vol.h5m");
510 }
511
512 // Destroy DM, no longer needed.
513 CHKERR DMDestroy(&dm);
514
515 // Destroy ghost vector
516 CHKERR VecDestroy(&global_error);
517 }
519
520 // finish work cleaning memory, getting statistics, etc.
522
523 return 0;
524}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
static char help[]
constexpr double a
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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 ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
@ F
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:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:576
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:525
FTensor::Index< 'i', SPACE_DIM > i
double D
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
static const bool debug
PetscErrorCode DMMoFEMGetSubRowIS(DM dm, IS *is)
get sub problem is
Definition DMMoFEM.cpp:330
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
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:118
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
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:261
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition Simple.hpp:396
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
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:355
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
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:393
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:736
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:389
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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.