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 {
 
 
   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 .
#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 degrees of freedom on entity.
default operator for TRI element
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 data field.
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.