test_cache_on_entities.cpp File Reference

#include <MoFEM.hpp>

Go to the source code of this file.

Classes
struct	MyStorage
struct	OpVolumeSet
	Operator set cache stored data, in this example, global indices, but it can be any structure. More...
struct	OpVolumeTest
	Test if cached data can be accessed, and check consistency of data. More...
struct	OpVolumeAssemble
struct	VolRule
	Set integration rule. More...

Typedefs
using	VolEle = VolumeElementForcesAndSourcesCore
using	VolOp = VolEle::UserDataOperator
using	FaceEle = FaceElementForcesAndSourcesCore
using	FaceOp = FaceEle::UserDataOperator

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

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

Typedef Documentation

FaceEle

using FaceEle = FaceElementForcesAndSourcesCore

Definition at line 24 of file test_cache_on_entities.cpp.

FaceOp

using FaceOp = FaceEle::UserDataOperator

Definition at line 25 of file test_cache_on_entities.cpp.

VolEle

using VolEle = VolumeElementForcesAndSourcesCore

Examples

EshelbianPlasticity.cpp.

Definition at line 22 of file test_cache_on_entities.cpp.

VolOp

using VolOp = VolEle::UserDataOperator

Examples

EshelbianPlasticity.cpp, build_large_problem.cpp, and lorentz_force.cpp.

Definition at line 23 of file test_cache_on_entities.cpp.

Function Documentation

main()

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

Definition at line 167 of file test_cache_on_entities.cpp.

                                 {
 
  // initialize petsc
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    // Create MoAB database
    moab::Core moab_core;
    moab::Interface &moab = moab_core;
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core mofem_core(moab);
    MoFEM::Interface &m_field = mofem_core;
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Simple interface
    Simple *simple_interface;
    CHKERR m_field.getInterface(simple_interface);
    {
      // get options from command line
      CHKERR simple_interface->getOptions();
      // load mesh file
      CHKERR simple_interface->loadFile();
      // add fields
      CHKERR simple_interface->addDomainField("FIELD", H1,
                                              AINSWORTH_LEGENDRE_BASE, 3);
      // set fields order
      CHKERR simple_interface->setFieldOrder("FIELD", 4);
      // setup problem
      CHKERR simple_interface->setUp();
 
      // get dm
      auto dm = simple_interface->getDM();
      // create elements
      auto domain_fe = boost::make_shared<VolEle>(m_field);
      // set integration rule
      domain_fe->getRuleHook = VolRule();
 
      auto get_skin = [&]() {
        Range ents;
        CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, ents, true);
        Skinner skin(&m_field.get_moab());
        Range skin_ents;
        CHKERR skin.find_skin(0, ents, false, skin_ents);
        // filter not owned entities, those are not on boundary
        ParallelComm *pcomm =
            ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
        Range proc_skin;
        CHKERR pcomm->filter_pstatus(skin_ents,
                                     PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                     PSTATUS_NOT, -1, &proc_skin);
        Range adj;
        for (auto d : {0, 1, 2})
          CHKERR m_field.get_moab().get_adjacencies(proc_skin, d, false, adj,
                                                    moab::Interface::UNION);
        // proc_skin.merge(adj);
        return proc_skin;
      };
 
      auto get_mark_skin_dofs = [&](Range &&skin) {
        auto problem_manager = m_field.getInterface<ProblemsManager>();
        auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
        problem_manager->markDofs(simple_interface->getProblemName(), ROW, skin,
                                  *marker_ptr);
        return marker_ptr;
      };
 
      SmartPetscObj<Vec> v;
      CHKERR DMCreateGlobalVector_MoFEM(dm, v);
      auto mark_dofs = get_mark_skin_dofs(get_skin());
 
      // set operator to the volume elements
      domain_fe->getOpPtrVector().push_back(new OpVolumeSet("FIELD"));
      domain_fe->getOpPtrVector().push_back(new OpVolumeTest("FIELD"));
 
      domain_fe->getOpPtrVector().push_back(
          new OpSetBc("FIELD", true, mark_dofs));
      domain_fe->getOpPtrVector().push_back(new OpVolumeAssemble("FIELD", v));
      domain_fe->getOpPtrVector().push_back(new OpUnSetBc("FIELD"));
 
      // make integration in volume (here real calculations starts)
      CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                      domain_fe);
      OpVolumeSet::entsIndices.clear();
 
      MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
      CHKERR VecAssemblyBegin(v);
      CHKERR VecAssemblyEnd(v);
      CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
 
      CHKERR DMoFEMMeshToLocalVector(simple_interface->getDM(), v,
                                     INSERT_VALUES, SCATTER_REVERSE);
 
      // CHKERR VecView(v,PETSC_VIEWER_STDOUT_WORLD);
 
      const MoFEM::Problem *problem_ptr;
      CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
 
      double *array;
      CHKERR VecGetArray(v, &array);
      for (auto dof : *(problem_ptr->getNumeredRowDofsPtr())) {
        auto loc_idx = dof->getPetscLocalDofIdx();
 
        if (dof->getFieldData() != array[loc_idx])
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "Data inconsistency");
 
        if ((*mark_dofs)[loc_idx]) {
          if (array[loc_idx] != 0)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Dof assembled, but it should be not");
        } else {
          if (array[loc_idx] == 0)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Dof not assembled, but it should be");
        }
      }
 
      CHKERR VecRestoreArray(v, &array);
    }
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

help

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

static

Definition at line 20 of file test_cache_on_entities.cpp.