20static char help[] =
"...\n\n";
36 double operator()(
const double x,
const double y,
const double z)
const {
37 return 1 + x * x + y * y + z * z * z;
46 const double z)
const {
68 double operator()(
const double x,
const double y,
const double z)
const {
73struct OpS :
public FaceElementForcesAndSourcesCore::UserDataOperator {
75 OpS(
const bool beta = 1)
108 if (row_side == col_side && row_type == col_type) {
147 double area = getArea() *
bEta;
149 auto t_w = getFTensor0IntegrationWeight();
155 const double alpha = t_w * area;
160 for (
int rr = 0; rr !=
nbRows; rr++) {
164 for (
int cc = 0; cc !=
nbCols; cc++) {
166 a += alpha * t_row_base * t_col_base;
188 const int *row_indices = &*row_data.
getIndices().data().begin();
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;
195 &*
locMat.data().begin(), ADD_VALUES);
201 &*
locMat.data().begin(), ADD_VALUES);
207int main(
int argc,
char *argv[]) {
212 moab::Core moab_core;
213 moab::Interface &moab = moab_core;
219 PetscBool flg_test = PETSC_FALSE;
220 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Poisson's problem options",
226 CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
227 &flg_test, PETSC_NULLPTR);
247 boost::shared_ptr<ForcesAndSourcesCore>
249 boost::shared_ptr<ForcesAndSourcesCore>
251 boost::shared_ptr<ForcesAndSourcesCore>
253 boost::shared_ptr<ForcesAndSourcesCore>
255 boost::shared_ptr<ForcesAndSourcesCore>
257 boost::shared_ptr<PoissonExample::PostProcFE>
259 boost::shared_ptr<ForcesAndSourcesCore> null;
260 boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
267 boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe,
false);
272 global_error, domain_error);
277 const double beta = 1;
278 boundary_penalty_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
281 boundary_penalty_lhs_fe->getOpPtrVector().push_back(
new OpS(beta));
282 boundary_rhs_fe->getOpPtrVector().push_back(
357 CHKERR DMCreateGlobalVector(dm, &
F);
362 DM dm_sub_KK, dm_sub_LU;
363 ublas::matrix<Mat> nested_matrices(2, 2);
364 ublas::vector<IS> nested_is(2);
366 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_KK);
367 CHKERR DMSetType(dm_sub_KK,
"DMMOFEM");
369 CHKERR DMSetFromOptions(dm_sub_KK);
377 CHKERR DMSetUp(dm_sub_KK);
379 CHKERR DMCreateMatrix(dm_sub_KK, &nested_matrices(0, 0));
380 domain_lhs_fe->ksp_B = nested_matrices(0, 0);
383 boundary_penalty_lhs_fe->ksp_B = nested_matrices(0, 0);
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);
391 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_LU);
392 CHKERR DMSetType(dm_sub_LU,
"DMMOFEM");
394 CHKERR DMSetFromOptions(dm_sub_LU);
400 CHKERR DMSetUp(dm_sub_LU);
402 CHKERR DMCreateMatrix(dm_sub_LU, &nested_matrices(1, 0));
403 boundary_lhs_fe->ksp_B = nested_matrices(1, 0);
406 CHKERR MatAssemblyBegin(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
407 CHKERR MatAssemblyEnd(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
409 CHKERR MatTranspose(nested_matrices(1, 0), MAT_INITIAL_MATRIX,
410 &nested_matrices(0, 1));
411 CHKERR DMDestroy(&dm_sub_LU);
413 domain_rhs_fe->ksp_f =
F;
416 boundary_rhs_fe->ksp_f =
F;
423 nested_matrices(1, 1) = PETSC_NULLPTR;
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;
434 cerr <<
"UU" << endl;
435 MatView(nested_matrices(0, 0), PETSC_VIEWER_DRAW_WORLD);
437 cerr <<
"LU" << endl;
438 MatView(nested_matrices(1, 0), PETSC_VIEWER_DRAW_WORLD);
440 cerr <<
"UL" << endl;
441 MatView(nested_matrices(0, 1), PETSC_VIEWER_DRAW_WORLD);
445 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is[0], 2, &nested_is[0],
446 &nested_matrices(0, 0), &A);
450 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
451 CHKERR KSPSetFromOptions(solver);
453 CHKERR KSPSetOperators(solver, A, A);
455 CHKERR KSPGetPC(solver, &pc);
456 PetscBool is_pcfs = PETSC_FALSE;
457 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
459 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[0]);
460 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[1]);
463 "Only works with pre-conditioner PCFIELDSPLIT");
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));
498 if (flg_test == PETSC_TRUE) {
509 post_proc_volume->writeFile(
"out_vol.h5m");
516 CHKERR VecDestroy(&global_error);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
PetscErrorCode DMMoFEMGetSubRowIS(DM dm, IS *is)
get sub problem is
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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.
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.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
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.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
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.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
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.
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.