gel_jacobian_test.cpp

Go to the documentation of this file.

/** \file gel_jacobian_test.cpp
  \brief Atom test testing calculation of element residual vectors and tangent matrices
  \ingroup gel
*/
 
/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <PostProcOnRefMesh.hpp>
#include <Projection10NodeCoordsOnField.hpp>
 
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <adolc/adolc.h>
#include <Gels.hpp>
using namespace GelModule;
 
#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>
#include <fstream>
#include <iostream>
namespace bio = boost::iostreams;
using bio::tee_device;
using bio::stream;
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
 
  try {
 
  moab::Core mb_instance;
  moab::Interface& moab = mb_instance;
 
  {
    PetscBool flg = PETSC_TRUE;
    char mesh_file_name[255];
    ierr = PetscOptionsGetString(PETSC_NULL,PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
    if(flg != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
    }
    const char *option;
    option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
    rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
    ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
    if(pcomm == NULL) pcomm =  new ParallelComm(&moab,PETSC_COMM_WORLD);
  }
 
  MoFEM::Core core(moab);
  MoFEM::Interface& m_field = core;
  BitRefLevel bit_level0;
  bit_level0.set(0);
  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
 
  // Define fields and finite elements
  {
 
    // Set approximation fields
    {
      // Seed all mesh entities to MoFEM database, those entities can be potentially used as finite elements
      // or as entities which carry some approximation field.
 
      bool check_if_spatial_field_exist = m_field.check_field("SPATIAL_POSITION");
      ierr = m_field.add_field("SPATIAL_POSITION",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
      ierr = m_field.add_field("CHEMICAL_LOAD",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
      ierr = m_field.add_field("HAT_EPS",L2,AINSWORTH_LEGENDRE_BASE,6,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
      ierr = m_field.add_field("SPATIAL_POSITION_DOT",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
      ierr = m_field.add_field("CHEMICAL_LOAD_DOT",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
      ierr = m_field.add_field("HAT_EPS_DOT",L2,AINSWORTH_LEGENDRE_BASE,6,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
 
      //Add field H1 space rank 3 to approximate geometry using hierarchical basis
      //For 10 node tetrahedral, before use, geometry is projected on that field (see below)
      ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
 
      //meshset consisting all entities in mesh
      EntityHandle root_set = moab.get_root_set();
      //add entities to field (root_mesh, i.e. on all mesh entities fields are approx.)
      ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"SPATIAL_POSITION"); CHKERRQ(ierr);
      ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"CHEMICAL_LOAD"); CHKERRQ(ierr);
      ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"HAT_EPS"); CHKERRQ(ierr);
      ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"SPATIAL_POSITION_DOT"); CHKERRQ(ierr);
      ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"CHEMICAL_LOAD_DOT"); CHKERRQ(ierr);
      ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"HAT_EPS_DOT"); CHKERRQ(ierr);
 
      PetscBool flg;
      PetscInt order;
      ierr = PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
      if(flg != PETSC_TRUE) {
        order = 2;
      }
      if(order < 2) {
        //SETERRQ()
      }
 
      ierr = m_field.set_field_order(root_set,MBTET,"SPATIAL_POSITION",order); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBTRI,"SPATIAL_POSITION",order); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBEDGE,"SPATIAL_POSITION",order); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBVERTEX,"SPATIAL_POSITION",1); CHKERRQ(ierr);
 
      ierr = m_field.set_field_order(root_set,MBTET,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBTRI,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBEDGE,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBVERTEX,"SPATIAL_POSITION_DOT",1); CHKERRQ(ierr);
 
      ierr = m_field.set_field_order(root_set,MBTET,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBTRI,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBEDGE,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBVERTEX,"CHEMICAL_LOAD",1); CHKERRQ(ierr);
 
      ierr = m_field.set_field_order(root_set,MBTET,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBTRI,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBEDGE,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBVERTEX,"CHEMICAL_LOAD_DOT",1); CHKERRQ(ierr);
 
      ierr = m_field.set_field_order(root_set,MBTET,"HAT_EPS",order-1); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBTET,"HAT_EPS_DOT",order-1); CHKERRQ(ierr);
 
      //gemetry approximation is set to 2nd oreder
      ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
      ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
      ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
      ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
      ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
 
      ierr = m_field.build_fields(); CHKERRQ(ierr);
 
      // Sett geometry approximation and initial spatial positions
      // 10 node tets
      if (!check_if_spatial_field_exist) {
        Projection10NodeCoordsOnField ent_method_material(m_field, "MESH_NODE_POSITIONS");
        ierr = m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material); CHKERRQ(ierr);
        Projection10NodeCoordsOnField ent_method_spatial(m_field, "SPATIAL_POSITION");
        ierr = m_field.loop_dofs("SPATIAL_POSITION", ent_method_spatial); CHKERRQ(ierr);
      }
 
    }
 
    //Set finite elements
    {
      ierr = m_field.add_finite_element("GEL_FE",MF_ZERO); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_row("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_col("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_row("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_col("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_row("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_col("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
 
      ierr = m_field.modify_finite_element_add_field_data("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_data("GEL_FE","SPATIAL_POSITION_DOT"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_data("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_data("GEL_FE","CHEMICAL_LOAD_DOT"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_data("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_data("GEL_FE","HAT_EPS_DOT"); CHKERRQ(ierr);
      ierr = m_field.modify_finite_element_add_field_data("GEL_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
      EntityHandle root_set = moab.get_root_set();
      ierr = m_field.add_ents_to_finite_element_by_type(root_set,MBTET,"GEL_FE"); CHKERRQ(ierr);
 
      //build finite elemnts
      ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
      //build adjacencies
      ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
    }
 
  }
 
  // Create gel instance
  Gel gel(m_field);
  {
 
    map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
    gel.constitutiveEquationPtr = boost::shared_ptr<Gel::ConstitutiveEquation<adouble> >(
      new Gel::ConstitutiveEquation<adouble>(gel.blockMaterialData)
    );
 
    // Set material parameters
    material_data[0].gAlpha = 1;
    material_data[0].vAlpha = 0.3;
    material_data[0].gBeta = 1;
    material_data[0].vBeta = 0.3;
    material_data[0].gBetaHat = 1;
    material_data[0].vBetaHat = 0.3;
    material_data[0].oMega = 1;
    material_data[0].vIscosity = 1;
    material_data[0].pErmeability = 2.;
 
    int def_block = 0;
    Tag th_block_id;
    rval = moab.tag_get_handle(
      "BLOCK_ID",1,MB_TYPE_INTEGER,th_block_id,MB_TAG_CREAT|MB_TAG_SPARSE,&def_block
    ); CHKERRQ_MOAB(rval);
 
    Gel::CommonData &common_data = gel.commonData;
    common_data.spatialPositionName = "SPATIAL_POSITION";
    common_data.spatialPositionNameDot = "SPATIAL_POSITION_DOT";
    common_data.muName = "CHEMICAL_LOAD";
    common_data.muNameDot = "CHEMICAL_LOAD_DOT";
    common_data.strainHatName = "HAT_EPS";
    common_data.strainHatNameDot = "HAT_EPS_DOT";
 
    Gel::GelFE *fe_ptr[] = { &gel.feRhs, &gel.feLhs };
    for(int ss = 0;ss<2;ss++) {
      fe_ptr[ss]->getOpPtrVector().push_back(new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION",common_data,false,true));
      fe_ptr[ss]->getOpPtrVector().push_back(new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION_DOT",common_data,false,true));
      fe_ptr[ss]->getOpPtrVector().push_back(new Gel::OpGetDataAtGaussPts("CHEMICAL_LOAD",common_data,true,true));
      fe_ptr[ss]->getOpPtrVector().push_back(new Gel::OpGetDataAtGaussPts("HAT_EPS",common_data,true,false,MBTET));
      fe_ptr[ss]->getOpPtrVector().push_back(new Gel::OpGetDataAtGaussPts("HAT_EPS_DOT",common_data,true,false,MBTET));
    }
 
    // attach tags for each recorder
    vector<int> tags;
    tags.push_back(1);
    tags.push_back(2);
    tags.push_back(3);
    tags.push_back(4);
 
    // Right hand side operators
    gel.feRhs.getOpPtrVector().push_back(
      new Gel::OpJacobian("SPATIAL_POSITION", tags,gel.constitutiveEquationPtr,gel.commonData,true,false)
    );
    gel.feRhs.getOpPtrVector().push_back(
      new Gel::OpRhsStressTotal(gel.commonData)
    );
    gel.feRhs.getOpPtrVector().push_back(
      new Gel::OpRhsSolventFlux(gel.commonData)
    );
    gel.feRhs.getOpPtrVector().push_back(
      new Gel::OpRhsSolventConcetrationDot(gel.commonData)
    );
    gel.feRhs.getOpPtrVector().push_back(
      new Gel::OpRhsStrainHat(gel.commonData)
    );
    // Left hand side operators
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpJacobian("SPATIAL_POSITION",tags,gel.constitutiveEquationPtr,gel.commonData,false,true)
    );
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpLhsdxdx(gel.commonData)
    );
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpLhsdxdMu(gel.commonData)
    );
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpLhsdxdStrainHat(gel.commonData)
    );
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpLhsdStrainHatdStrainHat(gel.commonData)
    );
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpLhsdStrainHatdx(gel.commonData)
    );
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpLhsdMudMu(gel.commonData)
    );
    gel.feLhs.getOpPtrVector().push_back(
      new Gel::OpLhsdMudx(gel.commonData)
    );
 
  }
 
  // Create dm instance
  DM dm;
  DMType dm_name = "DMGEL";
  {
    ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
    ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
    ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
    ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
    ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
    //add elements to dm
    ierr = DMMoFEMAddElement(dm,"GEL_FE"); CHKERRQ(ierr);
    ierr = DMSetUp(dm); CHKERRQ(ierr);
  }
 
  // Make calculations
  Mat M;
  Vec F,U_t;
  {
    ierr = DMCreateGlobalVector_MoFEM(dm,&F); CHKERRQ(ierr);
    ierr = DMCreateGlobalVector_MoFEM(dm,&U_t); CHKERRQ(ierr);
    ierr = DMCreateMatrix_MoFEM(dm,&M); CHKERRQ(ierr);
    ierr = VecZeroEntries(F); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecZeroEntries(U_t); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(U_t,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(U_t,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = MatZeroEntries(M); CHKERRQ(ierr);
    gel.feRhs.ts_F = F; // Set right hand side vector manually
    gel.feRhs.ts_u_t = U_t;
    gel.feRhs.ts_ctx = TSMethod::CTX_TSSETIFUNCTION;
    ierr = DMoFEMLoopFiniteElements(dm,"GEL_FE",&gel.feRhs); CHKERRQ(ierr);
    gel.feLhs.ts_B = M; // Set matrix M
    gel.feLhs.ts_a = 1.0; // Set time step parameter
    gel.feLhs.ts_ctx = TSMethod::CTX_TSSETIJACOBIAN;
    ierr = DMoFEMLoopFiniteElements(dm,"GEL_FE",&gel.feLhs); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
    ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  }
 
  // See results
  {
 
    PetscViewer viewer;
    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"gel_jacobian_test.txt",&viewer); CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
    //MatView(M,PETSC_VIEWER_DRAW_WORLD);
    //std::string wait;
    //std::cin >> wait;
 
  }
 
  double mnorm;
  ierr = MatNorm(M,NORM_1,&mnorm); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"mnorm  = %9.8e\n",mnorm); CHKERRQ(ierr);
 
  if(fabs(mnorm-6.54443978e+01)>1e-6) {
    SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
  }
 
  // Clean and destroy
  {
    ierr = VecDestroy(&F); CHKERRQ(ierr);
    ierr = VecDestroy(&U_t); CHKERRQ(ierr);
    ierr = MatDestroy(&M); CHKERRQ(ierr);
    ierr = DMDestroy(&dm); CHKERRQ(ierr);
  }
  
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}