Gels.hpp

Go to the documentation of this file.

/** \file Gels.hpp
  \brief Implementation of Gel finite element
  \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/>. */
 
#ifndef __GEL_HPP__
#define __GEL_HPP__
 
#ifndef WITH_ADOL_C
  #error "MoFEM need to be compiled with ADOL-C"
#endif
 
namespace GelModule {
 
/** \brief Implementation of Gel constitutive model
\ingroup gel
 
Implementation follows constitutive equation from:
 
VISCOELASTICITY AND POROELASTICITY IN ELASTOMERIC GELS
Acta Mechanica Solida Sinica, Vol. 25, No. 5,
Yuhang Hu Zhigang Suo
 
Note1: Following implementation is tailored for large strain analysis,
however in current version the engineering strain is used, i.e. assuming small deformation.
Moreover calculation of solvent flux and other physical quantities is with
assumption that those values are nonlinear.
 
Note2: Basic implementation is with linear constitutive equations, however
adding nonlinearities to model should be simple, since all tangent matrices
are calculated using adloc-c.
 
System of equations of equations in the form is assembled:
\f[
\left[
\begin{array}{ccc}
\mathbf{K}_{xx} & \mathbf{K}_{x\mu} & \mathbf{K}_{x\hat{\varepsilon}} \\
\mathbf{K}_{\mu x} & \mathbf{K}_{\mu \mu} & \mathbf{K}_{\mu \hat{\varepsilon}} \\
\mathbf{K}_{\hat{\varepsilon} x} & \mathbf{K}_{\hat{\varepsilon} \mu} & \mathbf{K}_{\hat{\varepsilon} \hat{\varepsilon}}
\end{array}
\right]
\left[
\begin{array}{c}
\mathbf{q}_{x}\\
\mathbf{q}_{\mu}\\
\mathbf{q}_{\hat{\varepsilon}}
\end{array}
\right]
=
\left[
\begin{array}{l}
\mathbf{f}^\textrm{internal}_{x}\\
\mathbf{f}^V_{\mu}+\mathbf{f}^J_{\mu}\\
\mathbf{r}_{\hat{\varepsilon}}
\end{array}
\right]
\f]
 
*/
struct Gel {
 
  enum TagEvaluate {
    STRESSTOTAL,
    SOLVENTFLUX,
    SOLVENTRATE,
    RESIDUALSTRAINHAT
  };
 
  MoFEM::Interface &mFiled;
 
  /** \brief Gel material parameters
  \ingroup gel
  */
  struct BlockMaterialData {
 
    // FIXME: Need to have hood description for each parameter and units
    double vAlpha;        ///< Poisson ratio spring alpha
    double gAlpha;        ///< Shear modulus spring alpha
    double vBeta;         ///< Poisson ratio spring beta
    double gBeta;         ///< Shear modulus spring beta
    double vBetaHat;      ///< Poisson ratio dashpot
    double gBetaHat;      ///< Shear modulus dashpot
    double pErmeability;  ///< Permeability of Gel
    double vIscosity;     ///< Viscosity of Gel
    double oMega;         ///< Volume per solvent molecule
    double mU0;           ///< equilibrated chemical potential/load
 
    Range tEts;           ///< Test in the block
 
  };
  map<int,BlockMaterialData> blockMaterialData;
 
  /** \brief Constitutive model functions
  \ingroup gel
 
  \image html gel_spring_daspot_model.png "Gel model" width=700px
 
  */
  template<typename TYPE>
  struct ConstitutiveEquation {
 
    int iD;
    map<int,BlockMaterialData> &dAta;
 
    ConstitutiveEquation(map<int,BlockMaterialData> &data):
    dAta(data) {
    }
    virtual ~ConstitutiveEquation() {}
 
    // Input
 
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > F;            ///< Gradient of deformation
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > FDot;         ///< Rate of gradient of deformation
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainHat;    ///< Internal variable, strain in dashpot beta
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainHatDot; ///< Internal variable, strain in dashpot beta
    TYPE mU;    ///< Solvent concentration
    ublas::vector<TYPE,ublas::bounded_array<TYPE,3> > gradientMu;   ///< Gradient of solvent concentration
 
    // Internal use
 
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > C;              ///< Cauchy deformation
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > gradientU;      ///< Gradient of displacements
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainTotal;    ///< Total strain applied at integration point
 
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressAlpha;    ///< Stress generated by spring alpha
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressBeta;     ///< Stress generated by spring beta
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > strainHatFlux;  ///< Rate of dashpot (beta) strain
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressBetaHat;  ///< Stress as result of volume change due to solvent concentration
 
    TYPE traceStrainTotal;
    TYPE traceStrainHat;
    TYPE traceStressBeta;
    TYPE traceStrainTotalDot;
 
    // Output
 
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > stressTotal;  ///< Total stress
    ublas::vector<TYPE,ublas::bounded_array<TYPE,9> > solventFlux;                   ///< Solvent flux
    TYPE solventConcentrationDot;                                                       ///< Volume rate change
    ublas::matrix<TYPE,ublas::row_major,ublas::bounded_array<TYPE,9> > residualStrainHat;  ///< Residual for calculation epsilon hat
 
    virtual PetscErrorCode calculateCauchyDefromationTensor() {
      PetscFunctionBegin;
      C.resize(3,3);
      noalias(C) = prod(trans(F),F);
      PetscFunctionReturn(0);
    }
 
    /** \brief Calculate total strain
 
    */
    virtual PetscErrorCode calculateStrainTotal() {
      PetscFunctionBegin;
      gradientU.resize(3,3,false);
      noalias(gradientU) = F;
      for(int ii = 0;ii<3;ii++) {
        gradientU(ii,ii) -= 1;
      }
      strainTotal.resize(3,3,false);
      noalias(strainTotal) = gradientU + trans(gradientU);
      strainTotal *= 0.5;
      PetscFunctionReturn(0);
    }
 
    /** \berief Calculate trace of rate of gradient of deformation
    */
    virtual PetscErrorCode calculateTraceStrainTotalDot() {
      PetscFunctionBegin;
      traceStrainTotalDot = 0;
      for(int ii = 0;ii<3;ii++) {
        traceStrainTotalDot += FDot(ii,ii);
      }
      PetscFunctionReturn(0);
    }
 
    /** \brief Calculate stress in spring alpha
 
    \f[
    \sigma^\alpha_{ij} = 2G^\alpha\left(\varepsilon_{ij} + \frac{v^\alpha}{1-2v^\alpha}\varepsilon_{kk}\delta_{ij}\right)
    \f]
 
    Note: In general implementation for large strain this function should calculate
    Piola-Kirchhoff Stress I.
 
    */
    virtual PetscErrorCode calculateStressAlpha() {
      PetscFunctionBegin;
      traceStrainTotal = 0;
      for(int ii = 0;ii<3;ii++) {
        traceStrainTotal += strainTotal(ii,ii);
      }
      stressAlpha.resize(3,3,false);
      double a = 2.0*dAta[iD].gAlpha;
      double b = a*(dAta[iD].vAlpha/(1.0-2.0*dAta[iD].vAlpha));
      noalias(stressAlpha) = a*strainTotal;
      for(int ii = 0; ii<3; ii++){
        stressAlpha(ii,ii) += b*traceStrainTotal;
      }
      PetscFunctionReturn(0);
    }
 
    /** \brief Calculate stress in spring beta
 
    \f[
    \sigma^\beta_{ij} = 2G^\beta\left[
    (\varepsilon_{ij}-\hat{\varepsilon}_{ij})
    + \frac{v^\beta}{1-2v^\beta}(\varepsilon_{kk}-\hat{\varepsilon}_{kk})\delta_{ij}
    \right]
    \f]
 
    */
    virtual PetscErrorCode calculateStressBeta() {
      PetscFunctionBegin;
      traceStrainTotal = 0;
      traceStrainHat = 0;
      for(int ii = 0;ii<3;ii++) {
        traceStrainHat += strainHat(ii,ii);
        traceStrainTotal += strainTotal(ii,ii);
      }
      stressBeta.resize(3,3,false);
      double a = 2.0*dAta[iD].gBeta;
      double b = a*(dAta[iD].vBeta/(1.0-2.0*dAta[iD].vBeta));
      noalias(stressBeta) = a*(strainTotal-strainHat);
      for(int ii = 0;ii<3;ii++) {
        stressBeta(ii,ii) += b*(traceStrainTotal-traceStrainHat);
      }
      PetscFunctionReturn(0);
    }
 
    /** \brief Calculate rate of strain hat
 
    \f[
    \frac{\partial \hat{\varepsilon}_{ij}}{\partial t} =
    \frac{1}{2\hat{G}^\beta}
    \left(
    \sigma_{ij}^\beta - \frac{\hat{v}^\beta}{1+\hat{v}^\beta}\sigma_{kk}^\beta\delta_{ij}
    \right)
    \f]
 
    */
    virtual PetscErrorCode calculateStrainHatFlux() {
      PetscFunctionBegin;
      traceStressBeta = 0;
      for(int ii = 0;ii<3;ii++) {
        traceStressBeta = stressBeta(ii,ii);
      }
      strainHatFlux.resize(3,3,false);
      double a = (1.0/(2.0*dAta[iD].gBetaHat));
      noalias(strainHatFlux) = a*stressBeta;
      double b = a*(dAta[iD].vBetaHat/(1.0+dAta[iD].vBetaHat));
      for(int ii = 0;ii<3;ii++) {
        strainHatFlux(ii,ii) -= b*traceStressBeta;
      }
      PetscFunctionReturn(0);
    }
 
    /** \brief Calculate stress due to concentration of solvent molecules
 
    \f[
    \hat{\sigma}_{ij}^\beta = \frac{\Delta\mu}{\Omega}\delta_{ij}
    \f]
    note that \f$\Delta\mu = \mu-\mu_0\f$. In current version of code
    value of \f$\Delta\mu\f$ is approximated.
 
    */
    virtual PetscErrorCode calculateStressBetaHat() {
      PetscFunctionBegin;
      stressBetaHat.resize(3,3,false);
      stressBetaHat.clear();
      for(int ii=0; ii<3; ii++){
        stressBetaHat(ii,ii) = -(mU-dAta[iD].mU0)/dAta[iD].oMega;
      }
      PetscFunctionReturn(0);
    }
 
    // Functions calculating output variables
 
    virtual PetscErrorCode calculateStressTotal() {
      PetscFunctionBegin;
      stressTotal.resize(3,3,false);
      noalias(stressTotal) = stressAlpha;
      noalias(stressTotal) += stressBeta;
      noalias(stressTotal) += stressBetaHat;
      PetscFunctionReturn(0);
    }
 
    virtual PetscErrorCode calculateResidualStrainHat() {
      PetscFunctionBegin;
      residualStrainHat.resize(3,3,false);
      noalias(residualStrainHat) =  strainHatDot - strainHatFlux;
      PetscFunctionReturn(0);
    }
 
    /** \brief Calculate flux
 
    \f[
    J_k =
    -\left(
    \frac{\kappa}{\eta\Omega^2}
    \right)
    \mu_{,i}
    \f]
 
    */
    virtual PetscErrorCode calculateSolventFlux() {
      PetscFunctionBegin;
      double a = dAta[iD].pErmeability/(dAta[iD].vIscosity*dAta[iD].oMega/**dAta[iD].oMega*/);
      solventFlux=a*gradientMu;
      PetscFunctionReturn(0);
    }
 
    /** \brief Calculate solvent concentration rate
 
      FIXME: For simplicity as first approximation set volume change
      as trace of gradient total strain
 
    */
    virtual PetscErrorCode calculateSolventConcentrationDot() {
      PetscFunctionBegin;
      solventConcentrationDot = traceStrainTotalDot/*/dAta[iD].oMega*/;
      PetscFunctionReturn(0);
    }
 
  };
 
  boost::shared_ptr<Gel::ConstitutiveEquation<adouble> > constitutiveEquationPtr;
 
  /** \brief Common data for gel model
  \ingroup gel
  */
  struct CommonData {
 
    string spatialPositionName;
    string spatialPositionNameDot;
    string strainHatName;
    string strainHatNameDot;
    string muName;
    string muNameDot;
 
    map<string,vector<VectorDouble> > dataAtGaussPts;
    map<string,vector<MatrixDouble> > gradAtGaussPts;
 
    vector<MatrixDouble3by3> stressTotal;
    vector<VectorDouble3> solventFlux;
    vector<double> solventConcentrationDot;
    vector<VectorDouble9> residualStrainHat;
 
    vector<double*> jacRowPtr;
    vector<MatrixDouble> jacStressTotal;
    vector<MatrixDouble> jacSolventFlux;
    vector<VectorDouble> jacSolventDot;
    vector<MatrixDouble> jacStrainHat;
 
    bool recordOn;
    map<int,int> nbActiveVariables,nbActiveResults;
 
    CommonData():
    recordOn(true) {
    }
 
  };
  CommonData commonData;
 
  /// \brief  definition of volume element
  struct GelFE: public MoFEM::VolumeElementForcesAndSourcesCore {
 
    CommonData &commonData;
    int addToRule; ///< Takes into account HO geometry
 
    GelFE(MoFEM::Interface &m_field,CommonData &common_data):
    MoFEM::VolumeElementForcesAndSourcesCore(m_field),
    commonData(common_data),
    addToRule(1) {
    }
 
    int getRule(int order) {
      return 2*(order-1)+addToRule;
    }
 
    PetscErrorCode preProcess() {
      
      PetscFunctionBegin;
      ierr = MoFEM::VolumeElementForcesAndSourcesCore::preProcess(); CHKERRQ(ierr);
 
      if(ts_ctx == CTX_TSSETIFUNCTION) {
 
        ierr = mField.set_other_local_ghost_vector(
          problemPtr,commonData.spatialPositionName,commonData.spatialPositionNameDot,COL,ts_u_t,INSERT_VALUES,SCATTER_REVERSE
        ); CHKERRQ(ierr);
        ierr = mField.set_other_local_ghost_vector(
          problemPtr,commonData.muName,commonData.muNameDot,COL,ts_u_t,INSERT_VALUES,SCATTER_REVERSE
        ); CHKERRQ(ierr);
        ierr = mField.set_other_local_ghost_vector(
          problemPtr,commonData.strainHatName,commonData.strainHatNameDot,COL,ts_u_t,INSERT_VALUES,SCATTER_REVERSE
        ); CHKERRQ(ierr);
 
      }
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode postProcess() {
      
      PetscFunctionBegin;
 
      ierr = MoFEM::VolumeElementForcesAndSourcesCore::postProcess(); CHKERRQ(ierr);
 
      if(ts_ctx == CTX_TSSETIFUNCTION) {
        ierr = VecAssemblyBegin(ts_F); CHKERRQ(ierr);
        ierr = VecAssemblyEnd(ts_F); CHKERRQ(ierr);
      }
      if(ts_ctx == CTX_TSSETIJACOBIAN) {
        ierr = MatAssemblyBegin(ts_B,MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(ts_B,MAT_FLUSH_ASSEMBLY); CHKERRQ(ierr);
      }
 
      PetscFunctionReturn(0);
    }
 
  };
 
  GelFE feRhs,feLhs;
 
  Gel(MoFEM::Interface &m_field):
  mFiled(m_field),
  feRhs(m_field,commonData),
  feLhs(m_field,commonData) {
  }
 
  struct OpGetDataAtGaussPts: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    bool calcVal;
    bool calcGrad;
    EntityType zeroAtType;
 
    OpGetDataAtGaussPts(
      const string field_name,
      CommonData &common_data,
      bool calc_val,
      bool calc_grad,
      EntityType zero_at_type = MBVERTEX
    ):
    MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
      field_name,UserDataOperator::OPCOL
    ),
    commonData(common_data),
    calcVal(calc_val),
    calcGrad(calc_grad),
    zeroAtType(zero_at_type) {
    }
 
    /** \brief Operator field value
    *
    */
    PetscErrorCode doWork(
      int side,EntityType type,DataForcesAndSurcesCore::EntData &data
    ){
      PetscFunctionBegin;
 
      try {
 
        int nb_dofs = data.getFieldData().size();
        if(nb_dofs == 0) {
          PetscFunctionReturn(0);
        }
        int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
        int nb_gauss_pts = data.getN().size1();
 
        // Initialize
        if(calcVal) {
          commonData.dataAtGaussPts[rowFieldName].resize(nb_gauss_pts);
          for(int gg = 0;gg<nb_gauss_pts;gg++) {
            commonData.dataAtGaussPts[rowFieldName][gg].resize(rank,false);
          }
        }
        if(calcGrad) {
          commonData.gradAtGaussPts[rowFieldName].resize(nb_gauss_pts);
          for(int gg = 0;gg<nb_gauss_pts;gg++) {
            commonData.gradAtGaussPts[rowFieldName][gg].resize(rank,3,false);
          }
        }
 
        // Zero values
        if(type == zeroAtType) {
          for(int gg = 0;gg<nb_gauss_pts;gg++) {
            if(calcVal) {
              commonData.dataAtGaussPts[rowFieldName][gg].clear();
            }
            if(calcGrad) {
              commonData.gradAtGaussPts[rowFieldName][gg].clear();
            }
          }
        }
 
        VectorDouble &values = data.getFieldData();
 
        // Calculate values at integration points
        for(int gg = 0;gg<nb_gauss_pts;gg++) {
          VectorDouble N = data.getN(gg,nb_dofs/rank);
          MatrixDouble diffN = data.getDiffN(gg,nb_dofs/rank);
          for(int dd = 0;dd<nb_dofs/rank;dd++) {
            for(int rr1 = 0;rr1<rank;rr1++) {
              if(calcVal) {
                commonData.dataAtGaussPts[rowFieldName][gg][rr1] += N[dd]*values[rank*dd+rr1];
              }
              if(calcGrad) {
                for(int rr2 = 0;rr2<3;rr2++) {
                  commonData.gradAtGaussPts[rowFieldName][gg](rr1,rr2) += diffN(dd,rr2)*values[rank*dd+rr1];
                }
              }
            }
          }
        }
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
 
      PetscFunctionReturn(0);
    }
 
  };
 
  struct OpJacobian: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
    vector<int> tagS;
    boost::shared_ptr<Gel::ConstitutiveEquation<adouble> > cE;
    CommonData &commonData;
 
    bool calculateResidualBool;
    bool calculateJacobianBool;
    bool &recordOn;
    map<int,int> &nbActiveVariables;
    map<int,int> &nbActiveResults;
 
    OpJacobian(
      const string field_name,
      vector<int> tags,
      boost::shared_ptr<Gel::ConstitutiveEquation<adouble> > ce,
      CommonData &common_data,
      bool calculate_residual,
      bool calculate_jacobian
    ):
    MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
      field_name,UserDataOperator::OPROW
    ),
    tagS(tags),
    cE(ce),
    commonData(common_data),
    calculateResidualBool(calculate_residual),
    calculateJacobianBool(calculate_jacobian),
    recordOn(common_data.recordOn),
    nbActiveVariables(common_data.nbActiveVariables),
    nbActiveResults(common_data.nbActiveResults) {
    }
 
    int nbGaussPts;
    VectorDouble activeVariables;
 
    PetscErrorCode recordStressTotal() {
      
      PetscFunctionBegin;
 
      try {
        if(tagS[STRESSTOTAL]<0) {
          PetscFunctionReturn(0);
        }
        cE->F.resize(3,3,false);
        cE->strainHat.resize(3,3,false);
        MatrixDouble &F = (commonData.gradAtGaussPts[commonData.spatialPositionName])[0];
        VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[0];
        VectorDouble mu = (commonData.dataAtGaussPts[commonData.muName])[0];
        trace_on(tagS[STRESSTOTAL]);
        {
          // Activate gradient of defamation
          nbActiveVariables[tagS[STRESSTOTAL]] = 0;
          for(int dd1 = 0;dd1<3;dd1++) {
            for(int dd2 = 0;dd2<3;dd2++) {
              cE->F(dd1,dd2) <<= F(dd1,dd2);
              nbActiveVariables[tagS[STRESSTOTAL]]++;
            }
          }
          // Using vector notation to store hatStrain
          // xx,yy,zz,2xy,2yz,2zy
          cE->strainHat(0,0) <<= strain_hat[0];
          cE->strainHat(1,1) <<= strain_hat[1];
          cE->strainHat(2,2) <<= strain_hat[2];
          cE->strainHat(0,1) <<= strain_hat[3];
          cE->strainHat(1,2) <<= strain_hat[4];
          cE->strainHat(0,2) <<= strain_hat[5];
          nbActiveVariables[tagS[STRESSTOTAL]] += 6;
          cE->strainHat(1,0) = cE->strainHat(0,1);
          cE->strainHat(2,1) = cE->strainHat(1,2);
          cE->strainHat(2,0) = cE->strainHat(0,2);
          cE->mU <<= mu[0];
          nbActiveVariables[tagS[STRESSTOTAL]]++;
          // Do calculations
          ierr = cE->calculateStrainTotal(); CHKERRQ(ierr);
          ierr = cE->calculateStressAlpha(); CHKERRQ(ierr);
          ierr = cE->calculateStressBeta(); CHKERRQ(ierr);
          ierr = cE->calculateStressBetaHat(); CHKERRQ(ierr);
          // Finally calculate result
          ierr = cE->calculateStressTotal(); CHKERRQ(ierr);
          // Results
          nbActiveResults[tagS[STRESSTOTAL]] = 0;
          commonData.stressTotal.resize(nbGaussPts);
          commonData.stressTotal[0].resize(3,3,false);
          for(int d1 = 0;d1<3;d1++) {
            for(int d2 = 0;d2<3;d2++) {
              cE->stressTotal(d1,d2) >>= (commonData.stressTotal[0])(d1,d2);
              nbActiveResults[tagS[STRESSTOTAL]]++;
            }
          }
        }
        trace_off();
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode recordSolventFlux() {
      PetscFunctionBegin;
      try {
        if(tagS[SOLVENTFLUX]<0) {
          PetscFunctionReturn(0);
        }
        
        cE->gradientMu.resize(3,false);
        MatrixDouble &gradient_mu = (commonData.gradAtGaussPts[commonData.muName])[0];
        trace_on(tagS[SOLVENTFLUX]);
        {
          // Activate rate of gradient of defamation
          nbActiveVariables[tagS[SOLVENTFLUX]] = 0;
          for(int ii = 0;ii<3;ii++) {
            cE->gradientMu[ii] <<= gradient_mu(0,ii);
            nbActiveVariables[tagS[SOLVENTFLUX]]++;
          }
          ierr = cE->calculateSolventFlux(); CHKERRQ(ierr);
          nbActiveResults[tagS[SOLVENTFLUX]] = 0;
          commonData.solventFlux.resize(nbGaussPts);
          commonData.solventFlux[0].resize(3,false);
          for(int d1 = 0;d1<3;d1++) {
            cE->solventFlux[d1] >>= commonData.solventFlux[0][d1];
            nbActiveResults[tagS[SOLVENTFLUX]]++;
          }
        }
        trace_off();
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode recordSolventConcentrationDot() {
      PetscFunctionBegin;
 
      if(tagS[SOLVENTRATE]<0) {
        PetscFunctionReturn(0);
      }
 
      
 
      MatrixDouble &F = (commonData.gradAtGaussPts[commonData.spatialPositionName])[0];
      MatrixDouble &F_dot = (commonData.gradAtGaussPts[commonData.spatialPositionNameDot])[0];
 
      trace_on(tagS[SOLVENTRATE]);
      {
        // Activate rate of gradient of defamation
        nbActiveVariables[tagS[SOLVENTRATE]] = 0;
        cE->F.resize(3,3,false);
        for(int dd1 = 0;dd1<3;dd1++) {
          for(int dd2 = 0;dd2<3;dd2++) {
            cE->F(dd1,dd2) <<= F(dd1,dd2);
            nbActiveVariables[tagS[SOLVENTRATE]]++;
          }
        }
        cE->FDot.resize(3,3,false);
        for(int dd1 = 0;dd1<3;dd1++) {
          for(int dd2 = 0;dd2<3;dd2++) {
            cE->FDot(dd1,dd2) <<= F_dot(dd1,dd2);
            nbActiveVariables[tagS[SOLVENTRATE]]++;
          }
        }
        ierr = cE->calculateTraceStrainTotalDot(); CHKERRQ(ierr);
        ierr = cE->calculateSolventConcentrationDot(); CHKERRQ(ierr);
        nbActiveResults[tagS[SOLVENTRATE]] = 0;
        commonData.solventConcentrationDot.resize(nbGaussPts);
        cE->solventConcentrationDot >>= commonData.solventConcentrationDot[0];
        nbActiveResults[tagS[SOLVENTRATE]]++;
      }
      trace_off();
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode recordResidualStrainHat() {
      PetscFunctionBegin;
 
      if(tagS[RESIDUALSTRAINHAT]<0) {
        PetscFunctionReturn(0);
      }
 
      
 
      cE->F.resize(3,3,false);
      cE->strainHat.resize(3,3,false);
      cE->strainHatDot.resize(3,3,false);
      MatrixDouble &F = (commonData.gradAtGaussPts[commonData.spatialPositionName])[0];
      VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[0];
      VectorDouble &strain_hat_dot = (commonData.dataAtGaussPts[commonData.strainHatNameDot])[0];
      trace_on(tagS[RESIDUALSTRAINHAT]);
      {
        // Activate gradient of defamation
        nbActiveVariables[tagS[RESIDUALSTRAINHAT]] = 0;
        for(int dd1 = 0;dd1<3;dd1++) {
          for(int dd2 = 0;dd2<3;dd2++) {
            cE->F(dd1,dd2) <<= F(dd1,dd2);
            nbActiveVariables[tagS[RESIDUALSTRAINHAT]]++;
          }
        }
        // Using vector notation to store hatStrain
        // xx,yy,zz,2xy,2yz,2zy
        cE->strainHat(0,0) <<= strain_hat[0];
        cE->strainHat(1,1) <<= strain_hat[1];
        cE->strainHat(2,2) <<= strain_hat[2];
        cE->strainHat(0,1) <<= strain_hat[3];
        cE->strainHat(1,2) <<= strain_hat[4];
        cE->strainHat(0,2) <<= strain_hat[5];
        nbActiveVariables[tagS[RESIDUALSTRAINHAT]] += 6;
        cE->strainHat(1,0) = cE->strainHat(0,1);
        cE->strainHat(2,1) = cE->strainHat(1,2);
        cE->strainHat(2,0) = cE->strainHat(0,2);
        // Activate strain hat dot
        cE->strainHatDot(0,0) <<= strain_hat_dot[0];
        cE->strainHatDot(1,1) <<= strain_hat_dot[1];
        cE->strainHatDot(2,2) <<= strain_hat_dot[2];
        cE->strainHatDot(0,1) <<= strain_hat_dot[3];
        cE->strainHatDot(1,2) <<= strain_hat_dot[4];
        cE->strainHatDot(0,2) <<= strain_hat_dot[5];
        nbActiveVariables[tagS[RESIDUALSTRAINHAT]] += 6;
        cE->strainHatDot(1,0) = cE->strainHatDot(0,1);
        cE->strainHatDot(2,1) = cE->strainHatDot(1,2);
        cE->strainHatDot(2,0) = cE->strainHatDot(0,2);
        ierr = cE->calculateStrainTotal(); CHKERRQ(ierr);
        ierr = cE->calculateStressBeta(); CHKERRQ(ierr);
        ierr = cE->calculateStrainHatFlux(); CHKERRQ(ierr);
        ierr = cE->calculateResidualStrainHat(); CHKERRQ(ierr);
        nbActiveResults[tagS[RESIDUALSTRAINHAT]] = 0;
        commonData.residualStrainHat.resize(nbGaussPts);
        commonData.residualStrainHat[0].resize(6,false);
        cE->residualStrainHat(0,0) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
        cE->residualStrainHat(1,1) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
        cE->residualStrainHat(2,2) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
        cE->residualStrainHat(0,1) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
        cE->residualStrainHat(1,2) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
        cE->residualStrainHat(0,2) >>= commonData.residualStrainHat[0][nbActiveResults[tagS[RESIDUALSTRAINHAT]]++];
      }
      trace_off();
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode calculateFunction(TagEvaluate te,double *ptr) {
      PetscFunctionBegin;
 
      int r;
      //play recorder for values
      r = ::function(
        tagS[te],
        nbActiveResults[tagS[te]],
        nbActiveVariables[tagS[te]],
        &activeVariables[0],
        ptr
      );
      if(r<3) { // function is locally analytic
        SETERRQ1(
          PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,
          "ADOL-C function evaluation with error r = %d",r
        );
      }
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode calculateJacobian(TagEvaluate te) {
      PetscFunctionBegin;
      try {
        int r;
        r = jacobian(
          tagS[te],
          nbActiveResults[tagS[te]],
          nbActiveVariables[tagS[te]],
          &activeVariables[0],
          &(commonData.jacRowPtr[0])
        );
        if(r<3) {
          SETERRQ(
            PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error"
          );
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode calculateAtIntPtsStressTotal() {
      PetscFunctionBegin;
 
      try {
 
        if(tagS[STRESSTOTAL]<0) {
          PetscFunctionReturn(0);
        }
 
 
        
        activeVariables.resize(nbActiveVariables[tagS[STRESSTOTAL]],false);
 
        for(int gg = 0;gg<nbGaussPts;gg++) {
 
          MatrixDouble &F = (commonData.gradAtGaussPts[commonData.spatialPositionName])[gg];
          VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[gg];
          double mu = (commonData.dataAtGaussPts[commonData.muName])[gg][0];
          int nb_active_variables = 0;
          // Activate gradient of defamation
          for(int dd1 = 0;dd1<3;dd1++) {
            for(int dd2 = 0;dd2<3;dd2++) {
              activeVariables[nb_active_variables++] = F(dd1,dd2);
            }
          }
          // Using vector notation to store hatStrain
          // xx,yy,zz,2xy,2yz,2zy
          for(int ii = 0;ii<6;ii++) {
            activeVariables[nb_active_variables++] = strain_hat[ii];
          }
          // chemical load
          activeVariables[nb_active_variables++] = mu;
 
          if(nb_active_variables!=nbActiveVariables[tagS[STRESSTOTAL]]) {
            SETERRQ(
              PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,
              "Number of active variables does not much"
            );
          }
 
          if(calculateResidualBool) {
            if(gg == 0) {
              commonData.stressTotal.resize(nbGaussPts);
            }
            commonData.stressTotal[gg].resize(3,3,false);
            ierr = calculateFunction(STRESSTOTAL,&commonData.stressTotal[gg](0,0)); CHKERRQ(ierr);
          }
 
          if(calculateJacobianBool) {
            if(gg == 0) {
              commonData.jacStressTotal.resize(nbGaussPts);
              commonData.jacRowPtr.resize(nbActiveResults[tagS[STRESSTOTAL]]);
            }
            commonData.jacStressTotal[gg].resize(
              nbActiveResults[tagS[STRESSTOTAL]],
              nbActiveVariables[tagS[STRESSTOTAL]],
              false
            );
            for(int dd = 0;dd<nbActiveResults[tagS[STRESSTOTAL]];dd++) {
              commonData.jacRowPtr[dd] = &commonData.jacStressTotal[gg](dd,0);
            }
            ierr = calculateJacobian(STRESSTOTAL); CHKERRQ(ierr);
          }
 
        }
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode calculateAtIntPtsSolventFlux() {
      PetscFunctionBegin;
 
      if(tagS[SOLVENTFLUX]<0) {
        PetscFunctionReturn(0);
      }
 
      
 
      activeVariables.resize(nbActiveResults[tagS[SOLVENTFLUX]]);
 
      for(int gg = 0;gg<nbGaussPts;gg++) {
        MatrixDouble &gradient_mu = (commonData.gradAtGaussPts[commonData.muName])[gg];
        int nb_active_variables = 0;
        // Activate rate of gradient of defamation
        for(int ii = 0;ii<3;ii++) {
          activeVariables[nb_active_variables++] = gradient_mu(0,ii);
        }
        if(nb_active_variables!=nbActiveVariables[tagS[SOLVENTFLUX]]) {
          SETERRQ(
            PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Number of active variables does not much"
          );
        }
        if(calculateResidualBool) {
          commonData.solventFlux[gg].resize(3,false);
          ierr = calculateFunction(
            SOLVENTFLUX,
            &(commonData.solventFlux[gg][0])
          ); CHKERRQ(ierr);
        }
        if(calculateJacobianBool) {
          if(gg == 0) {
            commonData.jacSolventFlux.resize(nbGaussPts);
            commonData.jacRowPtr.resize(nbActiveResults[tagS[SOLVENTFLUX]]);
          }
          commonData.jacSolventFlux[gg].resize(
            nbActiveResults[tagS[SOLVENTFLUX]],
            nbActiveVariables[tagS[SOLVENTFLUX]],
            false
          );
          for(int dd = 0;dd<nbActiveResults[tagS[SOLVENTFLUX]];dd++) {
            commonData.jacRowPtr[dd] = &(commonData.jacSolventFlux[gg](dd,0));
          }
          ierr = calculateJacobian(SOLVENTFLUX); CHKERRQ(ierr);
        }
      }
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode calculateAtIntPtsSolventDot() {
      PetscFunctionBegin;
 
      if(tagS[SOLVENTRATE]<0) {
        PetscFunctionReturn(0);
      }
 
      
 
      activeVariables.resize(nbActiveVariables[tagS[SOLVENTRATE]],false);
 
      for(int gg = 0;gg<nbGaussPts;gg++) {
        MatrixDouble &F = (commonData.gradAtGaussPts[commonData.spatialPositionName])[gg];
        MatrixDouble &F_dot = (commonData.gradAtGaussPts[commonData.spatialPositionNameDot])[gg];
        int nb_active_variables = 0;
        // Activate gradient of defamation
        for(int dd1 = 0;dd1<3;dd1++) {
          for(int dd2 = 0;dd2<3;dd2++) {
            activeVariables[nb_active_variables++] = F(dd1,dd2);
          }
        }
        for(int dd1 = 0;dd1<3;dd1++) {
          for(int dd2 = 0;dd2<3;dd2++) {
            activeVariables[nb_active_variables++] = F_dot(dd1,dd2);
          }
        }
        if(nb_active_variables!=nbActiveVariables[tagS[SOLVENTRATE]]) {
          SETERRQ(
            PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Number of active variables does not much"
          );
        }
        if(calculateResidualBool) {
          ierr = calculateFunction(SOLVENTRATE,&commonData.solventConcentrationDot[gg]); CHKERRQ(ierr);
        }
        if(calculateJacobianBool) {
          if(gg == 0) {
            commonData.jacSolventDot.resize(nbGaussPts);
            commonData.jacRowPtr.resize(nbActiveResults[tagS[SOLVENTRATE]]);
          }
          commonData.jacSolventDot[gg].resize(nbActiveVariables[tagS[SOLVENTRATE]]);
          commonData.jacRowPtr[0] = &commonData.jacSolventDot[gg][0];
          ierr = calculateJacobian(SOLVENTRATE); CHKERRQ(ierr);
        }
      }
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode calculateAtIntPtrsResidualStrainHat() {
      PetscFunctionBegin;
 
      try {
        if(tagS[RESIDUALSTRAINHAT]<0) {
          PetscFunctionReturn(0);
        }
        
        activeVariables.resize(nbActiveVariables[tagS[RESIDUALSTRAINHAT]]);
        for(int gg = 0;gg<nbGaussPts;gg++) {
          MatrixDouble &F = (commonData.gradAtGaussPts[commonData.spatialPositionName])[gg];
          VectorDouble &strain_hat = (commonData.dataAtGaussPts[commonData.strainHatName])[gg];
          VectorDouble &strain_hat_dot = (commonData.dataAtGaussPts[commonData.strainHatNameDot])[gg];
          int nb_active_variables = 0;
          // Activate gradient of defamation
          for(int dd1 = 0;dd1<3;dd1++) {
            for(int dd2 = 0;dd2<3;dd2++) {
              activeVariables[nb_active_variables++] = F(dd1,dd2);
            }
          }
          // Using vector notation to store hatStrain
          // xx,yy,zz,2xy,2yz,2zy
          for(int ii = 0;ii<6;ii++) {
            activeVariables[nb_active_variables++] = strain_hat[ii];
          }
          // Using vector notation to store hatStrainDot
          for(int ii = 0;ii<6;ii++) {
            activeVariables[nb_active_variables++] = strain_hat_dot[ii];
          }
          if(nb_active_variables!=nbActiveVariables[tagS[RESIDUALSTRAINHAT]]) {
            SETERRQ(
              PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Number of active variables does not much"
            );
          }
          if(calculateResidualBool) {
            if(gg == 0) {
              commonData.residualStrainHat.resize(nbGaussPts);
            }
            commonData.residualStrainHat[gg].resize(6,false);
            ierr = calculateFunction(
              RESIDUALSTRAINHAT,&commonData.residualStrainHat[gg][0]
            ); CHKERRQ(ierr);
          }
          if(calculateJacobianBool) {
            if(gg == 0) {
              commonData.jacStrainHat.resize(nbGaussPts);
              commonData.jacRowPtr.resize(nbActiveResults[tagS[RESIDUALSTRAINHAT]]);
            }
            commonData.jacStrainHat[gg].resize(
              nbActiveResults[tagS[RESIDUALSTRAINHAT]],
              nbActiveVariables[tagS[RESIDUALSTRAINHAT]],
              false
            );
            for(int dd = 0;dd<nbActiveResults[tagS[RESIDUALSTRAINHAT]];dd++) {
              commonData.jacRowPtr[dd] = &commonData.jacStrainHat[gg](dd,0);
            }
            ierr = calculateJacobian(RESIDUALSTRAINHAT); CHKERRQ(ierr);
          }
        }
 
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
 
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode doWork(
      int row_side,EntityType row_type,DataForcesAndSurcesCore::EntData &row_data
    ) {
      PetscFunctionBegin;
      
      
 
      if(row_type != MBVERTEX) PetscFunctionReturn(0);
      nbGaussPts = row_data.getN().size1();
 
      Tag th_block_id;
      EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
      rval = getVolumeFE()->
      mField.get_moab().tag_get_handle("BLOCK_ID",th_block_id); CHKERRQ_MOAB(rval);
      rval = getVolumeFE()->
      mField.get_moab().tag_get_data(th_block_id,&ent,1,&(cE->iD)); CHKERRQ_MOAB(rval);
 
      try {
 
        if(recordOn) {
          ierr = recordStressTotal(); CHKERRQ(ierr);
          ierr = recordSolventFlux(); CHKERRQ(ierr);
          ierr = recordSolventConcentrationDot(); CHKERRQ(ierr);
          ierr = recordResidualStrainHat(); CHKERRQ(ierr);
        }
 
        ierr = calculateAtIntPtsStressTotal(); CHKERRQ(ierr);
        ierr = calculateAtIntPtsSolventFlux(); CHKERRQ(ierr);
        ierr = calculateAtIntPtsSolventDot(); CHKERRQ(ierr);
        ierr = calculateAtIntPtrsResidualStrainHat(); CHKERRQ(ierr);
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
 
      PetscFunctionReturn(0);
    }
 
  };
 
  struct AssembleVector: MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    AssembleVector(string field_name):
    MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW) {
    }
    
    VectorDouble nF;
    PetscErrorCode aSemble(
      int row_side,EntityType row_type,DataForcesAndSurcesCore::EntData &row_data
    ) {
      PetscFunctionBegin;
      try {
        int nb_dofs = row_data.getIndices().size();
        int *indices_ptr = &row_data.getIndices()[0];
        ierr = VecSetValues(
          getFEMethod()->ts_F,nb_dofs,indices_ptr,&nF[0],ADD_VALUES
        ); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Assemble internal force vector
  \ingroup gel
 
  \f[
  (\mathbf{f}^\textrm{internal}_x)_i =
  \int_V
  \frac{\partial N_i}{\partial X_j} \sigma_{ij}
  \textrm{d}V
  \f]
 
  */
  struct OpRhsStressTotal: public AssembleVector {
    CommonData &commonData;
    OpRhsStressTotal(CommonData &common_data):
    AssembleVector(common_data.spatialPositionName),
    commonData(common_data) {
    }
    PetscErrorCode doWork(
      int row_side,EntityType row_type,DataForcesAndSurcesCore::EntData &row_data
    ) {
      PetscFunctionBegin;
      try {
        int nb_dofs = row_data.getIndices().size();
        if(!nb_dofs) {
          PetscFunctionReturn(0);
        }
        nF.resize(nb_dofs,false);
        nF.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_dofs/3);
          const MatrixDouble& stress = commonData.stressTotal[gg];
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          for(int dd = 0;dd<nb_dofs/3;dd++) {
            for(int rr = 0;rr<3;rr++) {
              for(int nn = 0;nn<3;nn++) {
                nF[3*dd+rr] += val*diffN(dd,nn)*stress(rr,nn);
              }
            }
          }
        }
        ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Calculate internal forces for solvent flux
  \ingroup gel
 
  \f[
  (\mathbf{f}^J_{\mu})_j =
  \int_V
  \frac{\partial N_j}{\partial X_k} J_k
  \textrm{d}V
  \f]
 
  */
  struct OpRhsSolventFlux: public AssembleVector {
    CommonData &commonData;
    OpRhsSolventFlux(CommonData &common_data):
    AssembleVector(common_data.muName),
    commonData(common_data) {
    }
    PetscErrorCode doWork(
      int row_side,EntityType row_type,DataForcesAndSurcesCore::EntData &row_data
    ) {
      PetscFunctionBegin;
      try {
        int nb_dofs = row_data.getIndices().size();
        if(!nb_dofs) {
          PetscFunctionReturn(0);
        }
        nF.resize(nb_dofs,false);
        nF.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_dofs);
          const VectorDouble& flux = commonData.solventFlux[gg];
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          for(int dd = 0;dd<nb_dofs;dd++) {
            for(int jj = 0;jj<3;jj++) {
              nF[dd] += val*diffN(dd,jj)*flux(jj);
            }
          }
        }
        ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Calculating right hand side
  \ingroup gel
 
  \f[
  (\mathbf{f}^V_\mu)_j =
  \int_V
  N_j \frac{\partial V}{\partial t}
  \textrm{d}V
  \f]
 
  */
  struct OpRhsSolventConcetrationDot: public AssembleVector {
    CommonData &commonData;
    OpRhsSolventConcetrationDot(CommonData &common_data):
    AssembleVector(common_data.muName),
    commonData(common_data) {
    }
    PetscErrorCode doWork(
      int row_side,EntityType row_type,DataForcesAndSurcesCore::EntData &row_data
    ) {
      PetscFunctionBegin;
      try {
        int nb_dofs = row_data.getIndices().size();
        if(!nb_dofs) {
          PetscFunctionReturn(0);
        }
        nF.resize(nb_dofs,false);
        nF.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          const VectorAdaptor &N = row_data.getN(gg,nb_dofs);
          double solvent_dot_dot = commonData.solventConcentrationDot[gg];
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          nF += val*N*solvent_dot_dot;
        }
        ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Residual strain hat
  \ingroup gel
 
  \f[
  (\mathbf{r}_{\hat{\varepsilon}})_j =
  \int_V
  N_j \left(
  \frac{\partial \hat{\varepsilon}}{\partial t}-
  f(\sigma^\beta)
  \right)
  \textrm{d}V
  \f]
 
  */
  struct OpRhsStrainHat: public AssembleVector {
    CommonData &commonData;
    OpRhsStrainHat(CommonData &common_data):
    AssembleVector(common_data.strainHatName),
    commonData(common_data) {
    }
    PetscErrorCode doWork(
      int row_side,EntityType row_type,DataForcesAndSurcesCore::EntData &row_data
    ) {
      PetscFunctionBegin;
      try {
        int nb_dofs = row_data.getIndices().size();
        if(!nb_dofs) {
          PetscFunctionReturn(0);
        }
        nF.resize(nb_dofs,false);
        nF.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          const VectorAdaptor &N = row_data.getN(gg,nb_dofs/6);
          VectorDouble9& residual_strain_hat = commonData.residualStrainHat[gg];
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          for(int dd = 0;dd<nb_dofs/6;dd++) {
            for(int rr = 0;rr<6;rr++) {
              nF[dd*6+rr] += val*N[dd]*residual_strain_hat[rr];
            }
          }
        }
        ierr = aSemble(row_side,row_type,row_data); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  struct AssembleMatrix: MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
    AssembleMatrix(string row_name,string col_name):
    MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
      row_name,col_name,UserDataOperator::OPROWCOL) {
    }
    
    MatrixDouble K,transK;
    PetscErrorCode aSemble(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      try {
        int nb_row = row_data.getIndices().size();
        int nb_col = col_data.getIndices().size();
        int *row_indices_ptr = &row_data.getIndices()[0];
        int *col_indices_ptr = &col_data.getIndices()[0];
        // cerr << K << endl;
        ierr = MatSetValues(
          getFEMethod()->ts_B,
          nb_row,row_indices_ptr,
          nb_col,col_indices_ptr,
          &K(0,0),ADD_VALUES
        ); CHKERRQ(ierr);
        if(sYmm) {
          // Assemble of diagonal terms
          if(row_side != col_side || row_type != col_type) {
            transK.resize(nb_col,nb_row,false);
            noalias(transK) = trans(K);
            ierr = MatSetValues(
              getFEMethod()->ts_B,
              nb_col,col_indices_ptr,
              nb_row,row_indices_ptr,
              &transK(0,0),ADD_VALUES
            ); CHKERRQ(ierr);
          }
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Assemble matrix \f$\mathbf{K}_{xx}\f$
  */
  struct OpLhsdxdx: public AssembleMatrix {
    CommonData &commonData;
    OpLhsdxdx(CommonData &common_data):
    AssembleMatrix(
      common_data.spatialPositionName,
      common_data.spatialPositionName
    ),
    commonData(common_data) {
    }
    MatrixDouble dStress_dx;
    PetscErrorCode get_dStress_dx(
      DataForcesAndSurcesCore::EntData &col_data,int gg
    ) {
      PetscFunctionBegin;
      try {
        int nb_col = col_data.getIndices().size();
        dStress_dx.resize(9,nb_col,false);
        dStress_dx.clear();
        const MatrixAdaptor diffN = col_data.getDiffN(gg,nb_col/3);
        MatrixDouble &jac_stress = commonData.jacStressTotal[gg];
        // FIXME: This can be implemented more efficiently. At this stage is
        // efficiency bottle neck.
        for(int dd = 0;dd<nb_col/3;dd++) {  // DoFs in column
          for(int jj = 0;jj<3;jj++) {       // cont. DoFs in column
            double a = diffN(dd,jj);
            for(int rr = 0;rr<3;rr++) {     // Loop over dsigma_ii/dX_rr
              for(int ii = 0;ii<9;ii++) {   // ii represents components of stress tensor
                dStress_dx(ii,3*dd+rr) += jac_stress(ii,3*rr+jj)*a;
              }
            }
          }
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      int nb_row = row_data.getIndices().size();
      int nb_col = col_data.getIndices().size();
      if(nb_row == 0) PetscFunctionReturn(0);
      if(nb_col == 0) PetscFunctionReturn(0);
      try {
        K.resize(nb_row,nb_col,false);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ierr = get_dStress_dx(col_data,gg); CHKERRQ(ierr);
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          //cerr << dStress_dx << endl;
          dStress_dx *= val;
          const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_row/3);
          { //integrate element stiffness matrix
            for(int dd1 = 0;dd1<nb_row/3;dd1++) {
              for(int rr1 = 0;rr1<3;rr1++) {
                for(int dd2 = 0;dd2<nb_col/3;dd2++) {
                  for(int rr2 = 0;rr2<3;rr2++) {
                    K(3*dd1+rr1,3*dd2+rr2) += (
                      diffN(dd1,0)*dStress_dx(3*rr1+0,3*dd2+rr2)+
                      diffN(dd1,1)*dStress_dx(3*rr1+1,3*dd2+rr2)+
                      diffN(dd1,2)*dStress_dx(3*rr1+2,3*dd2+rr2)
                    );
                  }
                }
              }
            }
          }
        }
        //cerr << "G " << getMoFEMFEPtr()->get_ref_ent() << endl << K << endl;
        ierr = aSemble(
          row_side,col_side,row_type,col_type,row_data,col_data
        ); CHKERRQ(ierr);
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Assemble matrix \f$\mathbf{K}_{x\mu}\f$
  */
  struct OpLhsdxdMu: public AssembleMatrix {
    CommonData &commonData;
    OpLhsdxdMu(CommonData &common_data):
    AssembleMatrix(
      common_data.spatialPositionName,common_data.muName
    ),
    commonData(common_data) {
      sYmm = false;
    }
    MatrixDouble dStress_dMu;
    PetscErrorCode get_dStress_dx(
      DataForcesAndSurcesCore::EntData &col_data,int gg
    ) {
      PetscFunctionBegin;
      try {
        int nb_col = col_data.getIndices().size();
        dStress_dMu.resize(9,nb_col,false);
        dStress_dMu.clear();
        const VectorAdaptor N = col_data.getN(gg);
        MatrixDouble &jac_stress = commonData.jacStressTotal[gg];
        for(int dd = 0;dd<nb_col;dd++) {
          double a = N[dd];
          for(int ii = 0;ii<9;ii++) {
            dStress_dMu(ii,dd) += jac_stress(ii,9+6+0)*a;
          }
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      int nb_row = row_data.getIndices().size();
      int nb_col = col_data.getIndices().size();
      if(nb_row == 0) PetscFunctionReturn(0);
      if(nb_col == 0) PetscFunctionReturn(0);
      try {
        K.resize(nb_row,nb_col,false);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ierr = get_dStress_dx(col_data,gg); CHKERRQ(ierr);
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          dStress_dMu *= val;
          const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_row/3);
          { //integrate element stiffness matrix
            for(int dd1 = 0;dd1<nb_row/3;dd1++) {
              for(int rr1 = 0;rr1<3;rr1++) {
                for(int dd2 = 0;dd2<nb_col;dd2++) {
                  K(3*dd1+rr1,dd2) +=
                  diffN(dd1,0)*dStress_dMu(3*rr1+0,dd2)+
                  diffN(dd1,1)*dStress_dMu(3*rr1+1,dd2)+
                  diffN(dd1,2)*dStress_dMu(3*rr1+2,dd2);
                }
              }
            }
          }
        }
        ierr = aSemble(
          row_side,col_side,row_type,col_type,row_data,col_data
        ); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Assemble matrix \f$\mathbf{K}_{x\hat{\varepsilon}}\f$
  */
  struct OpLhsdxdStrainHat: public AssembleMatrix {
    CommonData &commonData;
    OpLhsdxdStrainHat(CommonData &common_data):
    AssembleMatrix(
      common_data.spatialPositionName,common_data.strainHatName
    ),
    commonData(common_data) {
      sYmm = false;
    }
    MatrixDouble dStress_dStrainHat;
    PetscErrorCode get_dStress_dStrainHat(
      DataForcesAndSurcesCore::EntData &col_data,int gg
    ) {
      PetscFunctionBegin;
      try {
        int nb_col = col_data.getIndices().size();
        dStress_dStrainHat.resize(9,nb_col,false);
        dStress_dStrainHat.clear();
        const VectorAdaptor N = col_data.getN(gg);
        MatrixDouble &jac_stress = commonData.jacStressTotal[gg];
        /*cerr << N << endl;
        cerr << jac_stress << endl;
        cerr << dStress_dStrainHat << endl;
        cerr << nb_col << endl;*/
        for(int dd = 0;dd<nb_col/6;dd++) {  /// DoFS in column
          double a = N[dd];
          for(int ii = 0;ii<9;ii++) {   // ii for elements in stress matrix
            for(int rr = 0;rr<6;rr++) {
              dStress_dStrainHat(ii,6*dd+rr) += jac_stress(ii,9+rr)*a;
            }
          }
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      int nb_row = row_data.getIndices().size();
      int nb_col = col_data.getIndices().size();
      if(nb_row == 0) PetscFunctionReturn(0);
      if(nb_col == 0) PetscFunctionReturn(0);
      try {
        K.resize(nb_row,nb_col,false);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ierr = get_dStress_dStrainHat(col_data,gg); CHKERRQ(ierr);
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          dStress_dStrainHat *= val;
          const MatrixAdaptor &diffN = row_data.getDiffN(gg,nb_row/3);
          { //integrate element stiffness matrix
            for(int dd1 = 0;dd1<nb_row/3;dd1++) {
              for(int rr1 = 0;rr1<3;rr1++) {
                for(int dd2 = 0;dd2<nb_col/6;dd2++) {
                  for(int rr2 = 0;rr2<6;rr2++) {
                    K(3*dd1+rr1,6*dd2+rr2) +=
                    diffN(dd1,0)*dStress_dStrainHat(3*rr1+0,6*dd2+rr2)+
                    diffN(dd1,1)*dStress_dStrainHat(3*rr1+1,6*dd2+rr2)+
                    diffN(dd1,2)*dStress_dStrainHat(3*rr1+2,6*dd2+rr2);
                  }
                }
              }
            }
          }
        }
        ierr = aSemble(
          row_side,col_side,row_type,col_type,row_data,col_data
        ); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Assemble matrix \f$\mathbf{K}_{\hat{\varepsilon}\hat{\varepsilon}}\f$
  */
  struct OpLhsdStrainHatdStrainHat: public AssembleMatrix {
    CommonData &commonData;
    OpLhsdStrainHatdStrainHat(CommonData &common_data):
    AssembleMatrix(
      common_data.strainHatName,common_data.strainHatName
    ),
    commonData(common_data) {
    }
    MatrixDouble dStrainHat_dStrainHat;
    PetscErrorCode get_dStrainHat_dStrainHat(
      DataForcesAndSurcesCore::EntData &col_data,int gg
    ) {
      PetscFunctionBegin;
      try {
        int nb_col = col_data.getIndices().size();
        dStrainHat_dStrainHat.resize(6,nb_col,false);
        dStrainHat_dStrainHat.clear();
        const VectorAdaptor N = col_data.getN(gg);
        MatrixDouble &jac_res_strain_hat = commonData.jacStrainHat[gg];
        for(int dd = 0;dd<nb_col/6;dd++) {  /// DoFS in column
          double a = N[dd];
          for(int ii = 0;ii<6;ii++) {   // ii for elements in stress matrix
            for(int rr = 0;rr<6;rr++) {
              dStrainHat_dStrainHat(ii,6*dd+rr) += jac_res_strain_hat(ii,9+rr)*a;
              dStrainHat_dStrainHat(ii,6*dd+rr) += jac_res_strain_hat(ii,9+6+rr)*a*getFEMethod()->ts_a; // Time dependent element
            }
          }
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      int nb_row = row_data.getIndices().size();
      int nb_col = col_data.getIndices().size();
      if(nb_row == 0) PetscFunctionReturn(0);
      if(nb_col == 0) PetscFunctionReturn(0);
      try {
        K.resize(nb_row,nb_col,false);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ierr = get_dStrainHat_dStrainHat(col_data,gg); CHKERRQ(ierr);
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          dStrainHat_dStrainHat *= val;
          const VectorAdaptor &N = row_data.getN(gg);
          { //integrate element stiffness matrix
            for(int dd1 = 0;dd1<nb_row/6;dd1++) {
              for(int rr1 = 0;rr1<6;rr1++) {
                for(int dd2 = 0;dd2<nb_col;dd2++) {
                  K(6*dd1+rr1,dd2) += N[dd1]*dStrainHat_dStrainHat(rr1,dd2);
                }
              }
            }
          }
        }
        ierr = aSemble(
          row_side,col_side,row_type,col_type,row_data,col_data
        ); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Assemble matrix \f$\mathbf{K}_{\hat{\varepsilon}x}\f$
  */
  struct OpLhsdStrainHatdx: public AssembleMatrix {
    CommonData &commonData;
    OpLhsdStrainHatdx(CommonData &common_data):
    AssembleMatrix(
      common_data.strainHatName,common_data.spatialPositionName
    ),
    commonData(common_data) {
      sYmm = false;
    }
    MatrixDouble dStrainHat_dx;
    PetscErrorCode get_dStrainHat_dx(
      DataForcesAndSurcesCore::EntData &col_data,int gg
    ) {
      PetscFunctionBegin;
      try {
        int nb_col = col_data.getIndices().size();
        dStrainHat_dx.resize(6,nb_col,false);
        dStrainHat_dx.clear();
        const MatrixAdaptor diffN = col_data.getDiffN(gg,nb_col/3);
        MatrixDouble &jac_res_strain_hat = commonData.jacStrainHat[gg];
        //cerr << "a\n" << diffN << endl;
        //cerr << "b\n" << jac_res_strain_hat << endl;
        for(int dd = 0;dd<nb_col/3;dd++) {    // DoFs in column
          for(int jj = 0;jj<3;jj++) {      // cont. DoFs in column
            double a = diffN(dd,jj);
            for(int ii = 0;ii<6;ii++) {       // ii for elements in strain hat
              for(int rr2 = 0;rr2<3;rr2++) {
                dStrainHat_dx(ii,3*dd+rr2) += jac_res_strain_hat(ii,3*rr2+jj)*a;
              }
            }
          }
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      int nb_row = row_data.getIndices().size();
      int nb_col = col_data.getIndices().size();
      if(nb_row == 0) PetscFunctionReturn(0);
      if(nb_col == 0) PetscFunctionReturn(0);
      try {
        K.resize(nb_row,nb_col,false);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ierr = get_dStrainHat_dx(col_data,gg); CHKERRQ(ierr);
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          dStrainHat_dx *= val;
          //cerr << dStrainHat_dx << endl;
          const VectorAdaptor &N = row_data.getN(gg);
          { //integrate element stiffness matrix
            for(int dd1 = 0;dd1<nb_row/6;dd1++) {
              for(int rr1 = 0;rr1<6;rr1++) {
                for(int dd2 = 0;dd2<nb_col/3;dd2++) {
                  for(int rr2 = 0;rr2<3;rr2++) {
                    K(6*dd1+rr1,3*dd2+rr2) += N[dd1]*dStrainHat_dx(rr1,3*dd2+rr2);
                  }
                }
              }
            }
          }
        }
        ierr = aSemble(
          row_side,col_side,row_type,col_type,row_data,col_data
        ); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /** \brief Assemble matrix \f$\mathbf{K}_{\mu \mu}\f$
  */
  struct OpLhsdMudMu: public AssembleMatrix {
    CommonData &commonData;
    OpLhsdMudMu(CommonData &common_data):
    AssembleMatrix(
      common_data.muName,common_data.muName
    ),
    commonData(common_data) {
    }
    MatrixDouble dSolventFlux_dmu;
    PetscErrorCode get_dSolventFlux_dmu(
      DataForcesAndSurcesCore::EntData &col_data,int gg
    ) {
      PetscFunctionBegin;
      try {
        int nb_col = col_data.getIndices().size();
        dSolventFlux_dmu.resize(3,nb_col,false);
        dSolventFlux_dmu.clear();
        const MatrixAdaptor diffN = col_data.getDiffN(gg);
        MatrixDouble &jac_solvent_flux = commonData.jacSolventFlux[gg];
        for(int dd = 0;dd<nb_col;dd++) {  // DoFs in column
          for(int jj = 0;jj<3;jj++) {
            double a = diffN(dd,jj);
            for(int rr2 = 0;rr2<3;rr2++) {
              dSolventFlux_dmu(rr2,dd) += jac_solvent_flux(rr2,jj)*a;
            }
          }
        }
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      int nb_row = row_data.getIndices().size();
      int nb_col = col_data.getIndices().size();
      if(nb_row == 0) PetscFunctionReturn(0);
      if(nb_col == 0) PetscFunctionReturn(0);
      try {
        K.resize(nb_row,nb_col,false);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ierr = get_dSolventFlux_dmu(col_data,gg); CHKERRQ(ierr);
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          dSolventFlux_dmu *= val;
          //cerr << dSolventFlux_dmu << endl;
          const MatrixAdaptor &diffN = row_data.getDiffN(gg);
          { //integrate element stiffness matrix
            for(int dd1 = 0;dd1<nb_row;dd1++) {
              for(int dd2 = 0;dd2<nb_col;dd2++) {
                for(int jj = 0;jj<3;jj++) {
                  K(dd1,dd2) += diffN(dd1,jj)*dSolventFlux_dmu(jj,dd2);
                }
              }
            }
          }
        }
        ierr = aSemble(
          row_side,col_side,row_type,col_type,row_data,col_data
        ); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
 
  /** \brief Assemble matrix \f$\mathbf{K}_{\mu x}\f$
  */
  struct OpLhsdMudx: public AssembleMatrix {
    CommonData &commonData;
    OpLhsdMudx(CommonData &common_data):
    AssembleMatrix(
      common_data.muName,common_data.spatialPositionName
    ),
    commonData(common_data) {
      sYmm = false;
    }
    VectorDouble dSolventDot_dx;
    PetscErrorCode get_dSolventDot_dx(
      DataForcesAndSurcesCore::EntData &col_data,int gg
    ) {
      PetscFunctionBegin;
      try {
        int nb_col = col_data.getIndices().size();
        dSolventDot_dx.resize(nb_col,false);
        dSolventDot_dx.clear();
        const MatrixAdaptor diffN = col_data.getDiffN(gg);
        VectorDouble &jac_solvent_rate = commonData.jacSolventDot[gg];
        //cerr << dSolventDot_dx << endl;
        //cerr << jac_solvent_rate << endl;
        for(int dd = 0;dd<nb_col/3;dd++) {    // DoFs in column
          for(int jj = 0;jj<3;jj++) {
            double a = diffN(dd,jj);
            for(int rr2 = 0;rr2<3;rr2++) {
              dSolventDot_dx[3*dd+rr2] += jac_solvent_rate(0+3*rr2+jj)*a;
              dSolventDot_dx[3*dd+rr2] += jac_solvent_rate(9+3*rr2+jj)*a*getFEMethod()->ts_a;
            }
          }
        }
        //cerr << dSolventDot_dx << endl;
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
      int nb_row = row_data.getIndices().size();
      int nb_col = col_data.getIndices().size();
      if(nb_row == 0) PetscFunctionReturn(0);
      if(nb_col == 0) PetscFunctionReturn(0);
      try {
        K.resize(nb_row,nb_col,false);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ierr = get_dSolventDot_dx(col_data,gg); CHKERRQ(ierr);
          double val = getVolume()*getGaussPts()(3,gg);
          if(getHoGaussPtsDetJac().size()>0) {
            val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
          }
          dSolventDot_dx *= val;
          //cerr << dSolventDot_dx << endl;
          const VectorAdaptor &N = row_data.getN(gg);
          { //integrate element stiffness matrix
            for(int dd1 = 0;dd1<nb_row;dd1++) {
              for(int dd2 = 0;dd2<nb_col;dd2++) {
                K(dd1,dd2) += N[dd1]*dSolventDot_dx[dd2];
              }
            }
          }
        }
        ierr = aSemble(
          row_side,col_side,row_type,col_type,row_data,col_data
        ); CHKERRQ(ierr);
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
  };
 
  /**
   * \brief Used to post proc stresses, energy, source term.
   *
   * calculation of principal stresses. for e.g.
   *
   */
  struct OpPostProcGel: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    moab::Interface &postProcMesh;
    std::vector<EntityHandle> &mapGaussPts;
 
    OpPostProcGel(
      const string field_name,
      moab::Interface &post_proc_mesh,
      std::vector<EntityHandle> &map_gauss_pts,
      CommonData &common_data
    ):
    MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
    commonData(common_data),
    postProcMesh(post_proc_mesh),
    mapGaussPts(map_gauss_pts) {
      // Set that this operator is executed only for vertices on element
      doVertices = true;
      doEdges = false;
      doQuads = false;
      doTris = false;
      doTets = false;
      doPrisms = false;
    }
 
    PetscErrorCode doWork(
      int side,EntityType type,DataForcesAndSurcesCore::EntData &data
    ) {
      // 
      
      PetscFunctionBegin;
 
      if(type!=MBVERTEX) PetscFunctionReturn(9);
      double def_VAL[9];
      bzero(def_VAL,9*sizeof(double));
      Tag th_stress_total;
      rval = postProcMesh.tag_get_handle(
        "STRESS_TOTAL",9,MB_TYPE_DOUBLE,th_stress_total,MB_TAG_CREAT|MB_TAG_SPARSE,def_VAL
      ); CHKERRQ_MOAB(rval);
      Tag th_flux;
      rval = postProcMesh.tag_get_handle(
        "FLUX_CHEMICAL_LOAD",3,MB_TYPE_DOUBLE,th_flux,MB_TAG_CREAT|MB_TAG_SPARSE,def_VAL
      ); CHKERRQ_MOAB(rval);
 
      const int nb_gauss_pts = mapGaussPts.size();
      for(int gg = 0;gg!=nb_gauss_pts;gg++) {
        rval = postProcMesh.tag_set_data(th_stress_total,&mapGaussPts[gg],1,&(commonData.stressTotal[gg](0,0))); CHKERRQ_MOAB(rval);
        rval = postProcMesh.tag_set_data(th_flux,&mapGaussPts[gg],1,&(commonData.solventFlux[gg][0])); CHKERRQ_MOAB(rval);
      }
 
      PetscFunctionReturn(0);
 
    }
 
  };
 
 
  struct MonitorPostProc: public MoFEM::FEMethod {
 
    MoFEM::Interface &mField;
    string pRoblem;
    string fE;
    CommonData &commonData;
    boost::shared_ptr<Gel::ConstitutiveEquation<adouble> > &cE;
    vector<int> tAgs;
 
    PostProcVolumeOnRefinedMesh postProc;
 
    bool iNit;
    int pRT;
 
    MonitorPostProc(
      MoFEM::Interface &m_field,
      string problem,
      string fe,
      CommonData &common_data,
      boost::shared_ptr<Gel::ConstitutiveEquation<adouble> > &ce,
      vector<int> tags
    ):
    FEMethod(),
    mField(m_field),
    pRoblem(problem),
    fE(fe),
    commonData(common_data),
    cE(ce),
    tAgs(tags),
    postProc(m_field),
    iNit(false) {
      
      PetscBool flg = PETSC_TRUE;
      ierr = PetscOptionsGetInt(
        PETSC_NULL,PETSC_NULL,"-my_output_prt",&pRT,&flg
      ); CHKERRABORT(PETSC_COMM_WORLD,ierr);
      if(flg!=PETSC_TRUE) {
        pRT = 1;
      }
    }
 
    PetscErrorCode preProcess() {
      PetscFunctionBegin;
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode operator()() {
      PetscFunctionBegin;
      PetscFunctionReturn(0);
    }
 
    PetscErrorCode postProcess() {
      PetscFunctionBegin;
      
      if(!iNit) {
 
        ierr = postProc.generateReferenceElementMesh(); CHKERRQ(ierr);
        ierr = postProc.addFieldValuesPostProc(commonData.spatialPositionName); CHKERRQ(ierr);
        ierr = postProc.addFieldValuesGradientPostProc(commonData.spatialPositionName); CHKERRQ(ierr);
        ierr = postProc.addFieldValuesPostProc(commonData.muName); CHKERRQ(ierr);
        ierr = postProc.addFieldValuesPostProc(commonData.strainHatName); CHKERRQ(ierr);
        postProc.getOpPtrVector().push_back
        (new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION",commonData,false,true));
        postProc.getOpPtrVector().push_back
        (new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION_DOT",commonData,false,true));
        postProc.getOpPtrVector().push_back
        (new Gel::OpGetDataAtGaussPts("CHEMICAL_LOAD",commonData,true,true));
        postProc.getOpPtrVector().push_back
        (new Gel::OpGetDataAtGaussPts("HAT_EPS",commonData,true,false,MBTET));
        postProc.getOpPtrVector().push_back
        (new Gel::OpGetDataAtGaussPts("HAT_EPS_DOT",commonData,true,false,MBTET));
        postProc.getOpPtrVector().push_back
        (new Gel::OpJacobian("SPATIAL_POSITION",tAgs,cE,commonData,true,false));
        postProc.getOpPtrVector().push_back
        (new OpPostProcGel("SPATIAL_POSITION",postProc.postProcMesh,postProc.mapGaussPts,commonData));
 
 
        iNit = true;
      }
      int step;
      ierr = TSGetTimeStepNumber(ts,&step); CHKERRQ(ierr);
      
      if((step)%pRT==0) {
        ierr = mField.loop_finite_elements(pRoblem,fE,postProc); CHKERRQ(ierr);
        ostringstream sss;
        sss << "out_" << step << ".h5m";
        rval = postProc.postProcMesh.write_file(sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"); CHKERRQ_MOAB(rval);
      }
      PetscFunctionReturn(0);
    }
 
  };
 
};
 
}
 
#endif //__GEL_HPP__
 
/***************************************************************************//**
 * \defgroup gel Gel model
 * \ingroup user_modules
 ******************************************************************************/