analytical_poisson_field_split.cpp File Reference

#include <BasicFiniteElements.hpp>
#include <PoissonOperators.hpp>
#include <AuxPoissonFunctions.hpp>

Go to the source code of this file.

Classes
struct	ExactFunction
	Function. More...
struct	ExactFunctionGrad
	Exact gradient. More...
struct	ExactLaplacianFunction
	Laplacian of function. More...
struct	OpS

Functions
int	main (int argc, char *argv[])

Variables
static char	help [] = "...\n\n"
static const bool	debug = false

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

< Volume element for the matrix

< Boundary element for the matrix

< Volume element to assemble vector

< Volume element to assemble vector

< Volume element evaluate error

< Volume element to Post-process results

< Null element do nothing

Definition at line 207 of file analytical_poisson_field_split.cpp.

                                 {
 
  // Initialize PETSc
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
  // Create MoAB database
  moab::Core moab_core;              // create database
  moab::Interface &moab = moab_core; // create interface to database
 
  try {
 
    // Get command line options
    int order = 3;                    // default approximation order
    PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
                             "none");
    // Set approximation order
    CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
                           PETSC_NULL);
    // Set testing (used by CTest)
    CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
                            &flg_test, PETSC_NULL);
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core mofem_core(moab);           // create database
    MoFEM::Interface &m_field = mofem_core; // create interface to database
    // Register DM Manager
    CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
 
    // Create vector to store approximation global error
    Vec global_error;
    CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
 
    // First we crate elements, implementation of elements is problem
    // independent, we do not know yet what fields are present in the problem,
    // or even we do not decided yet what approximation base or spaces we are
    // going to use. Implementation of element is free from those constrains and
    // can be used in different context.
 
    // Elements used by KSP & DM to assemble system of equations
    boost::shared_ptr<ForcesAndSourcesCore>
        domain_lhs_fe; ///< Volume element for the matrix
    boost::shared_ptr<ForcesAndSourcesCore>
        boundary_lhs_fe; ///< Boundary element for the matrix
    boost::shared_ptr<ForcesAndSourcesCore>
        domain_rhs_fe; ///< Volume element to assemble vector
    boost::shared_ptr<ForcesAndSourcesCore>
        boundary_rhs_fe; ///< Volume element to assemble vector
    boost::shared_ptr<ForcesAndSourcesCore>
        domain_error; ///< Volume element evaluate error
    boost::shared_ptr<PoissonExample::PostProcFE>
        post_proc_volume; ///< Volume element to Post-process results
    boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
    boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
    {
      // Add problem specific operators the generic finite elements to calculate
      // matrices and vectors.
      CHKERR PoissonExample::CreateFiniteElements(m_field)
          .createFEToAssembleMatrixAndVector(
              ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
              boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe, false);
      // Add problem specific operators the generic finite elements to calculate
      // error on elements and global error in H1 norm
      CHKERR PoissonExample::CreateFiniteElements(m_field)
          .createFEToEvaluateError(ExactFunction(), ExactFunctionGrad(),
                                   global_error, domain_error);
      // Post-process results
      CHKERR PoissonExample::CreateFiniteElements(m_field)
          .creatFEToPostProcessResults(post_proc_volume);
 
      const double beta = 1;
      boundary_penalty_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
          new FaceElementForcesAndSourcesCore(m_field));
      boundary_penalty_lhs_fe->getRuleHook = PoissonExample::FaceRule();
      boundary_penalty_lhs_fe->getOpPtrVector().push_back(new OpS(beta));
      boundary_rhs_fe->getOpPtrVector().push_back(
          new PoissonExample::Op_g(ExactFunction(), "U", beta));
    }
 
    // Get simple interface is simplified version enabling quick and
    // easy construction of problem.
    Simple *simple_interface;
    // Query interface and get pointer to Simple interface
    CHKERR m_field.getInterface(simple_interface);
 
    // Build problem with simple interface
    {
 
      // Get options for simple interface from command line
      CHKERR simple_interface->getOptions();
      // Load mesh file to database
      CHKERR simple_interface->loadFile();
 
      // Add field on domain and boundary. Field is declared by space and base
      // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
      // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
      // number of coefficients for dof.
      //
      // Simple interface assumes that there are four types of field; 1) defined
      // on body domain, 2) fields defined on body boundary, 3) skeleton field
      // defined on finite element skeleton. Finally data field is defined on
      // body domain. Data field is not solved but used for post-process or to
      // keep material parameters, etc. Here we using it to calculate
      // approximation error on elements.
 
      // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
      // used to construct base functions.
      CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
                                              1);
      // Add Lagrange multiplier field on body boundary
      CHKERR simple_interface->addBoundaryField("L", H1,
                                                AINSWORTH_LEGENDRE_BASE, 1);
      // Add error (data) field, we need only L2 norm. Later order is set to 0,
      // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
      // element. You can use more DOFs and collate moments of error to drive
      // anisotropic h/p-adaptivity, however this is beyond this example.
      CHKERR simple_interface->addDataField("ERROR", L2,
                                            AINSWORTH_LEGENDRE_BASE, 1);
 
      // Set fields order domain and boundary fields.
      CHKERR simple_interface->setFieldOrder("U",
                                             order); // to approximate function
      CHKERR simple_interface->setFieldOrder("L",
                                             order); // to Lagrange multipliers
      CHKERR simple_interface->setFieldOrder(
          "ERROR", 0); // approximation order for error
 
      // Setup problem. At that point database is constructed, i.e. fields,
      // finite elements and problem data structures with local and global
      // indexing.
      CHKERR simple_interface->setUp();
    }
 
    // Get access to PETSC-MoFEM DM manager.
    // or more derails see
    // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
    // Form that point internal MoFEM data structures are linked with PETSc
    // interface. If DM functions contains string MoFEM is is MoFEM specific DM
    // interface function, otherwise DM function part of standard PETSc
    // interface.
 
    DM dm;
    // Get dm
    CHKERR simple_interface->getDM(&dm);
 
    // Solve problem, only PETEc interface here.
    {
 
      // Create the right hand side vector and vector of unknowns
      Vec F, D;
      CHKERR DMCreateGlobalVector(dm, &F);
      // Create unknown vector by creating duplicate copy of F vector. only
      // structure is duplicated no values.
      CHKERR VecDuplicate(F, &D);
 
      DM dm_sub_KK, dm_sub_LU;
      ublas::matrix<Mat> nested_matrices(2, 2);
      ublas::vector<IS> nested_is(2);
 
      CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_KK);
      CHKERR DMSetType(dm_sub_KK, "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(dm_sub_KK, dm, "SUB_KK");
      CHKERR DMSetFromOptions(dm_sub_KK);
      CHKERR DMMoFEMSetSquareProblem(dm_sub_KK, PETSC_TRUE);
      CHKERR DMMoFEMAddSubFieldRow(dm_sub_KK, "U");
      CHKERR DMMoFEMAddSubFieldCol(dm_sub_KK, "U");
      CHKERR DMMoFEMAddElement(dm_sub_KK,
                               simple_interface->getDomainFEName().c_str());
      CHKERR DMMoFEMAddElement(dm_sub_KK,
                               simple_interface->getBoundaryFEName().c_str());
      CHKERR DMSetUp(dm_sub_KK);
      CHKERR DMMoFEMGetSubRowIS(dm_sub_KK, &nested_is[0]);
      CHKERR DMCreateMatrix(dm_sub_KK, &nested_matrices(0, 0));
      domain_lhs_fe->ksp_B = nested_matrices(0, 0);
      CHKERR DMoFEMLoopFiniteElements(
          dm_sub_KK, simple_interface->getDomainFEName(), domain_lhs_fe);
      boundary_penalty_lhs_fe->ksp_B = nested_matrices(0, 0);
      CHKERR DMoFEMLoopFiniteElements(dm_sub_KK,
                                      simple_interface->getBoundaryFEName(),
                                      boundary_penalty_lhs_fe);
      CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
      CHKERR DMDestroy(&dm_sub_KK);
      //
      CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_LU);
      CHKERR DMSetType(dm_sub_LU, "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(dm_sub_LU, dm, "SUB_LU");
      CHKERR DMSetFromOptions(dm_sub_LU);
      CHKERR DMMoFEMSetSquareProblem(dm_sub_LU, PETSC_FALSE);
      CHKERR DMMoFEMAddSubFieldRow(dm_sub_LU, "L");
      CHKERR DMMoFEMAddSubFieldCol(dm_sub_LU, "U");
      CHKERR DMMoFEMAddElement(dm_sub_LU,
                               simple_interface->getBoundaryFEName().c_str());
      CHKERR DMSetUp(dm_sub_LU);
      CHKERR DMMoFEMGetSubRowIS(dm_sub_LU, &nested_is[1]);
      CHKERR DMCreateMatrix(dm_sub_LU, &nested_matrices(1, 0));
      boundary_lhs_fe->ksp_B = nested_matrices(1, 0);
      CHKERR DMoFEMLoopFiniteElements(
          dm_sub_LU, simple_interface->getBoundaryFEName(), boundary_lhs_fe);
      CHKERR MatAssemblyBegin(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
      // CHKERR MatCreateTranspose(nested_matrices(1,0),&nested_matrices(0,1));
      CHKERR MatTranspose(nested_matrices(1, 0), MAT_INITIAL_MATRIX,
                          &nested_matrices(0, 1));
      CHKERR DMDestroy(&dm_sub_LU);
 
      domain_rhs_fe->ksp_f = F;
      CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                      domain_rhs_fe);
      boundary_rhs_fe->ksp_f = F;
      CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
                                      boundary_rhs_fe);
      CHKERR VecAssemblyBegin(F);
      CHKERR VecAssemblyEnd(F);
 
      Mat A;
      nested_matrices(1, 1) = PETSC_NULL;
 
      if (debug) {
        MatType type;
        MatGetType(nested_matrices(0, 0), &type);
        cerr << "K " << type << endl;
        MatGetType(nested_matrices(1, 0), &type);
        cerr << "C " << type << endl;
        MatGetType(nested_matrices(0, 1), &type);
        cerr << "CT " << type << endl;
        std::string wait;
        cerr << "UU" << endl;
        MatView(nested_matrices(0, 0), PETSC_VIEWER_DRAW_WORLD);
        std::cin >> wait;
        cerr << "LU" << endl;
        MatView(nested_matrices(1, 0), PETSC_VIEWER_DRAW_WORLD);
        std::cin >> wait;
        cerr << "UL" << endl;
        MatView(nested_matrices(0, 1), PETSC_VIEWER_DRAW_WORLD);
        std::cin >> wait;
      }
 
      CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is[0], 2, &nested_is[0],
                           &nested_matrices(0, 0), &A);
 
      // Create solver and link it to DM
      KSP solver;
      CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
      CHKERR KSPSetFromOptions(solver);
      // Set operators
      CHKERR KSPSetOperators(solver, A, A);
      PC pc;
      CHKERR KSPGetPC(solver, &pc);
      PetscBool is_pcfs = PETSC_FALSE;
      PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
      if (is_pcfs) {
        CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[0]);
        CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[1]);
      } else {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Only works with pre-conditioner PCFIELDSPLIT");
      }
      // Set-up solver, is type of solver and pre-conditioners
      CHKERR KSPSetUp(solver);
      // At solution process, KSP solver using DM creates matrices, Calculate
      // values of the left hand side and the right hand side vector. then
      // solves system of equations. Results are stored in vector D.
      CHKERR KSPSolve(solver, F, D);
 
      // Scatter solution on the mesh. Stores unknown vector on field on the
      // mesh.
      CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
      // Clean data. Solver and vector are not needed any more.
      CHKERR KSPDestroy(&solver);
      for (int i = 0; i != 2; i++) {
        CHKERR ISDestroy(&nested_is[i]);
        for (int j = 0; j != 2; j++) {
          CHKERR MatDestroy(&nested_matrices(i, j));
        }
      }
      CHKERR MatDestroy(&A);
      CHKERR VecDestroy(&D);
      CHKERR VecDestroy(&F);
    }
 
    // Calculate error
    {
      // Loop over all elements in mesh, and run users operators on each
      // element.
      CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                      domain_error);
      CHKERR PoissonExample::AuxFunctions(m_field).assembleGhostVector(
          global_error);
      CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
      if (flg_test == PETSC_TRUE) {
        CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
      }
    }
 
    {
      // Loop over all elements in the mesh and for each execute
      // post_proc_volume element and operators on it.
      CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                      post_proc_volume);
      // Write results
      post_proc_volume->writeFile("out_vol.h5m");
    }
 
    // Destroy DM, no longer needed.
    CHKERR DMDestroy(&dm);
 
    // Destroy ghost vector
    CHKERR VecDestroy(&global_error);
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

debug

const bool debug = false

static

Definition at line 21 of file analytical_poisson_field_split.cpp.

help

char help[] = "...\n\n"

static

Definition at line 20 of file analytical_poisson_field_split.cpp.