gel_analysis.cpp

Go to the documentation of this file.

/** \file gel_analysis.cpp
  \brief Reads cubit file and solves problem with gel material
  \ingroup gel
 
  \todo Current version is limited only to one material. If one like to have
  general problem for nonlinear elasticity should implement general time
  dependent problem. If inertia terms need to be considered, this material
  should be add to nonlinear dynamics problem.
 
  \todo Internal history state variables need to be statically condensed. It
  can  be done  by implementing static condensation on finite element level or
  by implementing pre-conditioner.
 
*/
 
/* 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>
#include <UserGelModel.hpp>
using namespace GelModule;
 
#include <boost/program_options.hpp>
using namespace std;
namespace po = boost::program_options;
 
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <SurfacePressure.hpp>
#include <NodalForce.hpp>
#include <EdgeForce.hpp>
 
// Elements for applying fluxes, convection or radiation
// are used to apply solvent concentration.
#include <ThermalElement.hpp>
 
//See file users_modules/elasticity/TimeForceScale.hpp
#include <TimeForceScale.hpp>
 
#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";
 
struct BlockData {
  int oRder;
  BlockData():
  oRder(-1) {
  }
};
 
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_gel_config,flg_file;
  char mesh_file_name[255];
  char gel_config_file[255];
  int order = 2;
  PetscBool is_partitioned = PETSC_FALSE;
  PetscBool is_atom_test = PETSC_FALSE;
 
  {
 
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Elastic Config","none"); CHKERRQ(ierr);
    ierr = PetscOptionsString(
      "-my_file",
      "mesh file name","",
      "mesh.h5m",
      mesh_file_name,
      255,
      &flg_file
    ); CHKERRQ(ierr);
    ierr = PetscOptionsInt(
      "-my_order",
      "default approximation order",
      "",
      order,
      &order,
      PETSC_NULL
    ); CHKERRQ(ierr);
 
    ierr = PetscOptionsBool(
      "-my_is_partitioned",
      "set if mesh is partitioned (this result that each process keeps only part of the mes","",
      PETSC_FALSE,&is_partitioned,PETSC_NULL
    ); CHKERRQ(ierr);
    ierr = PetscOptionsString(
      "-my_gel_config",
      "gel configuration file name","",
      "gel_config.in",gel_config_file,255,&flg_gel_config
    ); CHKERRQ(ierr);
 
    ierr = PetscOptionsBool(
      "-my_is_atom_test",
      "is used with testing, exit with error when diverged","",
      PETSC_FALSE,&is_atom_test,PETSC_NULL
    ); CHKERRQ(ierr);
 
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);
 
    //Reade parameters from line command
    if(flg_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
    }
 
    if(is_partitioned == PETSC_TRUE) {
      //Read mesh to MOAB
      const char *option;
      option =
      "PARALLEL=BCAST_DELETE;"
      "PARALLEL_RESOLVE_SHARED_ENTS;"
      "PARTITION=PARALLEL_PARTITION;";
      rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
    } else {
      const char *option;
      option = "";
      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;
 
  // Seed all mesh entities to MoFEM database, those entities can be potentially used as finite elements
  // or as entities which carry some approximation field.
  BitRefLevel bit_level0;
  bit_level0.set(0);
  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
 
  // Define fields and finite elements
  Gel gel(m_field);
  map<int,ThermalElement::FluxData> set_of_solvent_fluxes;
  {
 
    // Set approximation fields
    {
 
      // Add fields
      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);
 
      // Set approximation order. Solvent concentration has one order less than
      // order of of spatial position field. Tests need to be maid if that is
      // good enough to have stability for this type of problem. If not bubble
      // functions could be easily added by increasing approximate order for
      // volume.
 
      ierr = m_field.set_field_order(root_set,MBTET,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBTRI,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
      ierr = m_field.set_field_order(root_set,MBEDGE,"SPATIAL_POSITION",order,2); 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);
 
      try {
 
        map<int,BlockData> block_data;
        map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
 
        ifstream ini_file(gel_config_file);
        po::variables_map vm;
        po::options_description config_file_options;
 
        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);
 
        for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
 
          const string &name = it->getName();
          if(name.compare(0,3,"GEL")!=0) continue;
 
          int block_id = it->getMeshsetId();
 
          ostringstream str_order;
          str_order << "block_" << block_id << ".oRder";
          config_file_options.add_options()
          (str_order.str().c_str(),po::value<int>(
            &block_data[block_id].oRder
          )->default_value(order));
          Range block_ents;
          rval = moab.get_entities_by_handle(
            it->getMeshset(),block_ents,true
          ); CHKERRQ_MOAB(rval);
          material_data[block_id].tEts = block_ents.subset_by_type(MBTET);
          std::vector<int> block_id_vec(material_data[block_id].tEts.size(),it->getMeshsetId());
          // cerr << "AAA " << block_id_vec.size() << " " << material_data[block_id].tEts.size() << endl;
          rval = moab.tag_set_data(
            th_block_id,material_data[block_id].tEts,&*block_id_vec.begin()
          );
 
          ostringstream str_g_alpha;
          str_g_alpha << "block_" << block_id << ".gAlpha";
          config_file_options.add_options()
          (str_g_alpha.str().c_str(),po::value<double>(&material_data[block_id].gAlpha)->default_value(1));
          ostringstream str_v_alpha;
          str_v_alpha << "block_" << block_id << ".vAlpha";
          config_file_options.add_options()
          (str_v_alpha.str().c_str(),po::value<double>(&material_data[block_id].vAlpha)->default_value(0));
          ostringstream str_g_beta;
          str_g_beta << "block_" << block_id << ".gBeta";
          config_file_options.add_options()
          (str_g_beta.str().c_str(),po::value<double>(&material_data[block_id].gBeta)->default_value(1));
          ostringstream str_v_beta;
          str_v_beta << "block_" << block_id << ".vBeta";
          config_file_options.add_options()
          (str_v_beta.str().c_str(),po::value<double>(&material_data[block_id].vBeta)->default_value(0));
          ostringstream str_g_beta_hat;
          str_g_beta_hat << "block_" << block_id << ".gBetaHat";
          config_file_options.add_options()
          (str_g_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].gBetaHat)->default_value(1));
          ostringstream str_v_beta_hat;
          str_v_beta_hat << "block_" << block_id << ".vBetaHat";
          config_file_options.add_options()
          (str_v_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].vBetaHat)->default_value(0));
          ostringstream str_omega;
          str_omega << "block_" << block_id << ".oMega";
          config_file_options.add_options()
          (str_omega.str().c_str(),po::value<double>(&material_data[block_id].oMega)->default_value(1));
          ostringstream str_viscosity;
          str_viscosity << "block_" << block_id << ".vIscosity";
          config_file_options.add_options()
          (str_viscosity.str().c_str(),po::value<double>(&material_data[block_id].vIscosity)->default_value(1));
          ostringstream str_permeability;
          str_permeability << "block_" << block_id << ".pErmeability";
          config_file_options.add_options()
          (str_permeability.str().c_str(),po::value<double>(&material_data[block_id].pErmeability)->default_value(1));
          ostringstream str_mu0;
          str_mu0 << "block_" << block_id << ".mU0";
          config_file_options.add_options()
          (str_mu0.str().c_str(),po::value<double>(&material_data[block_id].mU0)->default_value(1));
 
        }
 
        po::parsed_options parsed = parse_config_file(ini_file,config_file_options,true);
        store(parsed,vm);
        po::notify(vm);
        for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
 
          const string &name = it->getName();
          if(name.compare(0,3,"GEL")!=0) continue;
 
          if(block_data[it->getMeshsetId()].oRder == -1) continue;
          if(block_data[it->getMeshsetId()].oRder == order) continue;
          PetscPrintf(PETSC_COMM_WORLD,"Set block %d order to %d\n",it->getMeshsetId(),block_data[it->getMeshsetId()].oRder);
          Range block_ents;
          rval = moab.get_entities_by_handle(it->getMeshset(),block_ents,true); CHKERRQ_MOAB(rval);
          Range ents_to_set_order;
          rval = moab.get_adjacencies(block_ents,3,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
          ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
          rval = moab.get_adjacencies(block_ents,2,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
          rval = moab.get_adjacencies(block_ents,1,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
          ierr = m_field.synchronise_entities(ents_to_set_order); CHKERRQ(ierr);
          ierr = m_field.set_field_order(ents_to_set_order,"SPATIAL_POSITION",block_data[it->getMeshsetId()].oRder); CHKERRQ(ierr);
          ierr = m_field.set_field_order(ents_to_set_order,"CHEMICAL_LOAD",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
          ierr = m_field.set_field_order(ents_to_set_order,"HAT_EPS",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
          ierr = m_field.set_field_order(ents_to_set_order,"SPATIAL_POSITION_DOT",block_data[it->getMeshsetId()].oRder); CHKERRQ(ierr);
          ierr = m_field.set_field_order(ents_to_set_order,"CHEMICAL_LOAD_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
          ierr = m_field.set_field_order(ents_to_set_order,"HAT_EPS_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
        }
        vector<string> additional_parameters;
        additional_parameters = collect_unrecognized(parsed.options,po::include_positional);
        for(vector<string>::iterator vit = additional_parameters.begin();
        vit!=additional_parameters.end();vit++) {
          ierr = PetscPrintf(PETSC_COMM_WORLD,"** WARNING Unrecognised option %s\n",vit->c_str()); CHKERRQ(ierr);
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
      }
 
      for(
        map<int,Gel::BlockMaterialData>::iterator mit = gel.blockMaterialData.begin();
        mit!=gel.blockMaterialData.end();mit++
      ) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,"** Block Id %d\n",mit->first); CHKERRQ(ierr);
        PetscPrintf(PETSC_COMM_WORLD,"vAlpha = %5.4g\n",mit->second.vAlpha);
        PetscPrintf(PETSC_COMM_WORLD,"gAlpha = %5.4g\n",mit->second.gAlpha);
        PetscPrintf(PETSC_COMM_WORLD,"vBeta = %5.4g\n",mit->second.vBeta);
        PetscPrintf(PETSC_COMM_WORLD,"gBeta = %5.4g\n",mit->second.gBeta);
        PetscPrintf(PETSC_COMM_WORLD,"vBetaHat = %5.4g\n",mit->second.vBetaHat);
        PetscPrintf(PETSC_COMM_WORLD,"gBetaHat = %5.4g\n",mit->second.gBetaHat);
        PetscPrintf(PETSC_COMM_WORLD,"pErmeability = %5.4g\n",mit->second.pErmeability);
        PetscPrintf(PETSC_COMM_WORLD,"vIscosity = %5.4g\n",mit->second.vIscosity);
        PetscPrintf(PETSC_COMM_WORLD,"oMega = %5.4g\n",mit->second.oMega);
        PetscPrintf(PETSC_COMM_WORLD,"mU0 = %5.4g\n",mit->second.mU0);
      }
 
      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 value to chemical load
        map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
        for(
          map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
          mit!=material_data.end();mit++
        ) {
          // FIXME mu0 at border in undefined, if two blackset have different mu0
          Range vertices;
          rval = moab.get_connectivity(mit->second.tEts,vertices,true); CHKERRQ_MOAB(rval);
          ierr = m_field.set_field(mit->second.mU0,MBVERTEX,vertices,"CHEMICAL_LOAD"); CHKERRQ(ierr);
        }
      }
 
    }
 
    //Set finite elements. The primary element is GEL_FE in addition elements
    //for applying tractions and fluxes are added.
    {
      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);
 
      map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
      for(
        map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
        mit!=material_data.end();
        mit++
      ) {
        ierr = m_field.add_ents_to_finite_element_by_type(mit->second.tEts,MBTET,"GEL_FE"); CHKERRQ(ierr);
      }
 
      // Add Neumann forces, i.e. on triangles, edges and nodes.
      ierr = MetaNeummanForces::addNeumannBCElements(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
      ierr = MetaNodalForces::addElement(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
      ierr = MetaEdgeForces::addElement(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
 
      // Add solvent flux element
      {
 
        ierr = m_field.add_finite_element("CHEMICAL_LOAD_FLUX_FE",MF_ZERO); CHKERRQ(ierr);
        ierr = m_field.modify_finite_element_add_field_row("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
        ierr = m_field.modify_finite_element_add_field_col("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
        ierr = m_field.modify_finite_element_add_field_data("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
        ierr = m_field.modify_finite_element_add_field_data("CHEMICAL_LOAD_FLUX_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
        // Assume that boundary conditions are set in block containing surface
        // triangle elements and block name is "CHEMICAL_LOAD_FLUX"
        for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
          if(it->getName().compare(0,18,"FLUX_CHEMICAL_LOAD") == 0) {
            vector<double> data;
            ierr = it->getAttributes(data); CHKERRQ(ierr);
            if(data.size()!=1) {
              SETERRQ(PETSC_COMM_SELF,1,"Data inconsistency");
            }
            // Here it set how block of for heat flux is set.  This is because
            // implementation from thermal element is used to enforce this
            // boundary condition.
            strcpy(set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.name,"HeatFlu");
            set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.flag1 = 1;
            set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.value1 = data[0];
            //cerr << set_of_solvent_fluxes[it->getMeshsetId()].dAta << endl;
            rval = m_field.get_moab().get_entities_by_type(
              it->meshset,MBTRI,set_of_solvent_fluxes[it->getMeshsetId()].tRis,true
            ); CHKERRQ_MOAB(rval);
            ierr = m_field.add_ents_to_finite_element_by_type(
              set_of_solvent_fluxes[it->getMeshsetId()].tRis,MBTRI,"CHEMICAL_LOAD_FLUX_FE"
            ); CHKERRQ(ierr);
          }
        }
      }
 
      //build finite elements
      ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
      //build adjacencies
      ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
    }
 
  }
 
  // attach tags for each recorder
  vector<int> tags;
  tags.push_back(Gel::STRESSTOTAL); // ADOL-C tag used to calculate total stress
  tags.push_back(Gel::SOLVENTFLUX);
  tags.push_back(Gel::SOLVENTRATE);
  tags.push_back(Gel::RESIDUALSTRAINHAT);
 
  // Create gel instance and set operators.
  {
 
    gel.constitutiveEquationPtr = boost::shared_ptr<UserGelConstitutiveEquation<adouble> >(
      new UserGelConstitutiveEquation<adouble>(gel.blockMaterialData)
    );
 
    // Set name of fields which has been choose to approximate spatial
    // displacements, solvent concentration and internal state variables.
    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";
 
    // Set operators to calculate field values at integration points, both for
    // left and right hand side elements.
    Gel::GelFE *fe_ptr[] = { &gel.feRhs, &gel.feLhs };
    for(int ss = 0;ss<2;ss++) {
      ierr = addHOOpsVol("MESH_NODE_POSITIONS", *fe_ptr[ss], true, true, false,
                         false);
      CHKERRQ(ierr);
      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)
      );
    }
 
 
    // 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 discrete manager 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 = DMMoFEMAddElement(dm,"FORCE_FE"); CHKERRQ(ierr);
    ierr = DMMoFEMAddElement(dm,"PRESSURE_FE"); CHKERRQ(ierr);
    ierr = DMMoFEMAddElement(dm,"CHEMICAL_LOAD_FLUX_FE"); CHKERRQ(ierr);
    ierr = DMSetUp(dm); CHKERRQ(ierr);
  }
 
  // Create matrices and vectors used for analysis
  Vec T,F;
  Mat A;
  {
    ierr = DMCreateGlobalVector_MoFEM(dm,&T); CHKERRQ(ierr);
    ierr = VecDuplicate(T,&F); CHKERRQ(ierr);
    ierr = DMCreateMatrix_MoFEM(dm,&A); CHKERRQ(ierr);
  }
 
  // Setting finite element methods for Dirichelt boundary conditions
  SpatialPositionsBCFEMethodPreAndPostProc spatial_position_bc(
    m_field,"SPATIAL_POSITION",A,T,F
  );
  spatial_position_bc.methodsOp.push_back(new TimeForceScale("-my_displacements_history",false));
  DirichletBCFromBlockSetFEMethodPreAndPostProc concentration_bc(
    m_field,"CHEMICAL_LOAD","CHEMICAL_LOAD",A,T,F
  );
  concentration_bc.methodsOp.push_back(new TimeForceScale("-my_chemical_load_history",false));
 
  // Setting finite element method for applying tractions
  boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
  boost::ptr_map<string,NodalForce> nodal_forces;
  boost::ptr_map<string,EdgeForce> edge_forces;
  {
    //forces and pressures on surface
    ierr = MetaNeummanForces::setMomentumFluxOperators(m_field,neumann_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
    //noadl forces
    ierr = MetaNodalForces::setOperators(m_field,nodal_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
    //edge forces
    ierr = MetaEdgeForces::setOperators(m_field,edge_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
    for(
      boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
      mit!=neumann_forces.end();mit++
    ) {
      mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
    }
    for(
      boost::ptr_map<string,NodalForce>::iterator mit = nodal_forces.begin();
      mit!=nodal_forces.end();mit++
    ) {
      mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
    }
    for(
      boost::ptr_map<string,EdgeForce>::iterator mit = edge_forces.begin();
      mit!=edge_forces.end();mit++
    ) {
      mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
    }
  }
 
  // Solvent surface element, flux, convection radiation
  // TODO: add radiation and convection
  ThermalElement::MyTriFE solvent_surface_fe(m_field);
  {
    map<int,ThermalElement::FluxData>::iterator sit = set_of_solvent_fluxes.begin();
    for(;sit!=set_of_solvent_fluxes.end();sit++) {
      // add flux operator
      solvent_surface_fe.getOpPtrVector().push_back(
        new ThermalElement::OpHeatFlux("CHEMICAL_LOAD",PETSC_NULL,sit->second,true)
      );
    }
  }
 
  // Add finite elements to Time Stepping Solver, using Discrete Manager Interface
  {
    //Rhs
    ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIFunction(dm,"GEL_FE",&gel.feRhs,NULL,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIFunction(dm,"CHEMICAL_LOAD_FLUX_FE",&solvent_surface_fe,NULL,NULL); CHKERRQ(ierr);
    {
      boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
      for(;mit!=neumann_forces.end();mit++) {
        ierr = DMMoFEMTSSetIFunction(dm,mit->first.c_str(),&mit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
      }
    }
    {
      boost::ptr_map<string,NodalForce>::iterator fit = nodal_forces.begin();
      for(;fit!=nodal_forces.end();fit++) {
        ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
      }
    }
    {
      boost::ptr_map<string,EdgeForce>::iterator fit = edge_forces.begin();
      for(;fit!=edge_forces.end();fit++) {
        ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
      }
    }
    ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
 
    //Lhs
    ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIJacobian(dm,"GEL_FE",&gel.feLhs,NULL,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
    ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
 
  }
 
 
  // Create Time Stepping solver
  TS ts;
  {
    ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
    ierr = TSSetType(ts,TSBEULER); CHKERRQ(ierr);
  }
 
  {
    ierr = DMoFEMMeshToLocalVector(dm,T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = DMoFEMPreProcessFiniteElements(dm,&spatial_position_bc); CHKERRQ(ierr);
    ierr = DMoFEMPreProcessFiniteElements(dm,&concentration_bc); CHKERRQ(ierr);
    ierr = DMoFEMMeshToGlobalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
  }
 
  // Solve problem
  {
 
    ierr = TSSetIFunction(ts,F,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
    ierr = TSSetIJacobian(ts,A,A,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
 
    //Monitor
    Gel::MonitorPostProc post_proc(
      m_field,"DMGEL","GEL_FE",gel.commonData,gel.constitutiveEquationPtr,tags
    );
    TsCtx *ts_ctx;
    DMMoFEMGetTsCtx(dm,&ts_ctx);
    {
      ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
      ierr = TSMonitorSet(ts,f_TSMonitorSet,ts_ctx,PETSC_NULL); CHKERRQ(ierr);
    }
 
    double ftime = 1;
    ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime); CHKERRQ(ierr);
    ierr = TSSetFromOptions(ts); CHKERRQ(ierr);
    ierr = TSSetDM(ts,dm); CHKERRQ(ierr);
    #if PETSC_VERSION_GE(3,7,0)
    ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr);
    #endif
    ierr = TSSolve(ts,T); CHKERRQ(ierr);
    ierr = TSGetTime(ts,&ftime); CHKERRQ(ierr);
    PetscInt steps,snesfails,rejects,nonlinits,linits;
    ierr = TSGetTimeStepNumber(ts,&steps); CHKERRQ(ierr);
    ierr = TSGetSNESFailures(ts,&snesfails); CHKERRQ(ierr);
    ierr = TSGetStepRejections(ts,&rejects); CHKERRQ(ierr);
    ierr = TSGetSNESIterations(ts,&nonlinits); CHKERRQ(ierr);
    ierr = TSGetKSPIterations(ts,&linits); CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,
      "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",
      steps,rejects,snesfails,ftime,nonlinits,linits
    );
  }
 
  if(is_atom_test) {
    double sum = 0;
    ierr = VecSum(T,&sum); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"sum  = %9.8e\n",sum); CHKERRQ(ierr);
    double fnorm;
    ierr = VecNorm(T,NORM_2,&fnorm); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"fnorm  = %9.8e\n",fnorm); CHKERRQ(ierr);
 
    if (fabs(sum - 6.82084809e+01) > 1e-6) {
      SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
    }
    if (fabs(fnorm - 3.66184607e+01) > 1e-6) {
      SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
    }
  }
 
  // Clean and destroy
  {
    ierr = DMDestroy(&dm); CHKERRQ(ierr);
    ierr = VecDestroy(&T); CHKERRQ(ierr);
    ierr = VecDestroy(&F); CHKERRQ(ierr);
    ierr = MatDestroy(&A); CHKERRQ(ierr);
    //ierr = TSDestroy(&ts); CHKERRQ(ierr);
  }
 
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}