ComplexConstArea.hpp

Go to the documentation of this file.

/* 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 __COMPLEX_CONST_AREA_HPP__
#define __COMPLEX_CONST_AREA_HPP__
 
#include <math.h>
#include <complex>
 
extern "C" {
 
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3],
                        double circumcenter[3], double *xi, double *eta,
                        double *zeta);
void tricircumcenter3d_tp(double a[3], double b[3], double c[3],
                          double circumcenter[3], double *xi, double *eta);
}
 
namespace ObosleteUsersModules {
 
/**
 *
 * dN/dX = (1/A) * [ Spin[dX/dksi]*dN/deta - Spin[dX/deta]*dN/dksi ]
 *
 */
struct ConstrainConstantAarea : public FEMethod {
 
  MoFEM::Interface &mField;
 
  moab::Interface &moab;
 
  Mat C; //<< constrains matrix
  Mat Q; //<< projection matrix (usually this is shell matrix)
  string lambdaFieldName;
  int verbose;
 
  ConstrainConstantAarea(MoFEM::Interface &m_field, Mat _C, Mat _Q,
                         string lambda_field_name, int _verbose = 0)
      : mField(m_field), moab(m_field.get_moab()), C(_C), Q(_Q),
        lambdaFieldName(lambda_field_name), verbose(_verbose) {
    // calculate face shape functions direvatives
    diffNTRI.resize(3, 2);
    ShapeDiffMBTRI(&*diffNTRI.data().begin());
    // shape functions Gauss integration weigths
    G_TRI_W = G_TRI_W1;
    // nodal material positions
    dofs_X.resize(9);
    // noal values of Lagrange multipliers
    lambda.resize(3);
    // dofs indices for Lagrnage multipliers
    lambdaDofsRowIdx.resize(3);
    // lambda_dofs_row_ents.resize(3);
    lambdaDofsColIdx.resize(3);
    // dofs indices for rows and columns
    dispDofsRowIdx.resize(9);
    dispDofsColIdx.resize(9);
    localDispDofsRowIdx.resize(9);
    // face node coordinates
    coords.resize(9);
  }
 
  virtual ~ConstrainConstantAarea() = default;
 
  // elem data
  MatrixDouble diffNTRI;
  const double *G_TRI_W;
  vector<DofIdx> DirichletBC;
  VectorInt dispDofsColIdx, dispDofsRowIdx;
  VectorInt localDispDofsRowIdx;
  VectorInt lambdaDofsRowIdx, lambdaDofsColIdx;
  ublas::vector<double, ublas::bounded_array<double, 9>> coords;
  ublas::vector<double, ublas::bounded_array<double, 9>> dofs_X;
  ublas::vector<double, ublas::bounded_array<double, 3>> lambda;
  // vector<EntityHandle> lambda_dofs_row_ents;
  EntityHandle face;
 
  template <typename DOFS>
  inline auto loIt(DOFS dofs, const FieldBitNumber bit,
                   const EntityHandle ent) {
    return dofs->lower_bound(DofEntity::getLoFieldEntityUId(bit, ent));
  };
 
  template <typename DOFS>
  inline auto hiIt(DOFS dofs, const FieldBitNumber bit,
                   const EntityHandle ent) {
    return dofs->upper_bound(DofEntity::getHiFieldEntityUId(bit, ent));
  };
 
  /**
   * \brief get face data, indices and coors and nodal values
   *
   * \param is_that_C_otherwise_dC
   */
  PetscErrorCode getData(bool is_that_C_otherwise_dC, bool trans) {
    MoFEMFunctionBegin;
 
    face = numeredEntFiniteElementPtr->getEnt();
    fill(lambdaDofsRowIdx.begin(), lambdaDofsRowIdx.end(), -1);
    fill(lambdaDofsColIdx.begin(), lambdaDofsColIdx.end(), -1);
    fill(dispDofsRowIdx.begin(), dispDofsRowIdx.end(), -1);
    fill(dispDofsColIdx.begin(), dispDofsColIdx.end(), -1);
    fill(localDispDofsRowIdx.begin(), localDispDofsRowIdx.end(), -1);
    fill(lambda.begin(), lambda.end(), 0);
    // fill(lambda_dofs_row_ents.begin(),lambda_dofs_row_ents.end(),no_handle);
    const EntityHandle *conn_face;
    int num_nodes;
    CHKERR moab.get_connectivity(face, conn_face, num_nodes, true);
    if (num_nodes != 3)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "face should have three nodes");
    CHKERR moab.get_coords(conn_face, num_nodes, &*coords.data().begin());
 
    const auto field_bit_number = mField.get_field_bit_number(lambdaFieldName);
    const auto pos_bit_number =
        mField.get_field_bit_number("MESH_NODE_POSITIONS");
 
    auto row_dofs = getRowDofsPtr();
 
    for (int nn = 0; nn < num_nodes; nn++) {
 
      if (is_that_C_otherwise_dC) {
 
        auto dit = loIt(row_dofs, field_bit_number, conn_face[nn]);
        auto hi_dit = hiIt(row_dofs, field_bit_number, conn_face[nn]);
 
        // it is C
        // get rows which are Lagrange multipliers
        if (std::distance(dit, hi_dit) > 0) {
          if (std::distance(dit, hi_dit) != 1) {
            SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                     "data inconsistency, number of dof on node for < %s > "
                     "should be 1",
                     lambdaFieldName.c_str());
          }
          if (dit->get()->getPetscLocalDofIdx() < 0) {
            SETERRQ(
                PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "data inconsistency, negative index of local dofs on element");
          }
          lambda[nn] = dit->get()->getFieldData();
          lambdaDofsRowIdx[nn] = dit->get()->getPetscGlobalDofIdx();
          // lambda_dofs_row_ents[nn] = dit->getEnt();
        }
      }
      if ((!is_that_C_otherwise_dC) || (trans)) {
 
        auto dit = loIt(row_dofs, field_bit_number, conn_face[nn]);
        auto hi_dit = hiIt(row_dofs, field_bit_number, conn_face[nn]);
 
        // it is dC
        // get rows which are material nodal positions
        if (std::distance(dit, hi_dit) > 0) {
          if (std::distance(dit, hi_dit) != 1)
            SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                     "data inconsistency, number of dof on node for < %s > "
                     "should be 1",
                     lambdaFieldName.c_str());
          lambda[nn] = dit->get()->getFieldData();
 
          auto diit = loIt(row_dofs, pos_bit_number, conn_face[nn]);
          auto hi_diit = hiIt(row_dofs, pos_bit_number, conn_face[nn]);
 
          if (std::distance(diit, hi_diit) != 3)
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "data inconsistency, number of dof on node for "
                    "MESH_NODE_POSITIONS should be 3");
          for (; diit != hi_diit; diit++) {
            if (diit->get()->getPetscLocalDofIdx() < 0)
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                      "data inconsistency, negative index of local dofs on "
                      "element");
            assert(nn * 3 + diit->get()->getDofCoeffIdx() < 9);
            dispDofsRowIdx[nn * 3 + diit->get()->getDofCoeffIdx()] =
                diit->get()->getPetscGlobalDofIdx();
            localDispDofsRowIdx[nn * 3 + diit->get()->getDofCoeffIdx()] =
                diit->get()->getPetscLocalDofIdx();
          }
        }
      }
 
    
      auto col_dofs = getColDofsPtr();
 
      // get columns which are material nodal positions
      auto dit = loIt(col_dofs, pos_bit_number, conn_face[nn]);
      auto hi_dit = hiIt(col_dofs, pos_bit_number, conn_face[nn]);
      if (std::distance(dit, hi_dit) != 3) {
        SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "data inconsistency, number of dof on node for "
                 "MESH_NODE_POSITIONS should be 3, nb dofs = %d",
                 std::distance(dit, hi_dit));
      }
      for (; dit != hi_dit; dit++) {
        if (dit->get()->getPetscLocalDofIdx() < 0) {
          SETERRQ(
              PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "data inconsistency, negative index of local dofs on element");
        }
        dofs_X[nn * 3 + dit->get()->getDofCoeffIdx()] =
            dit->get()->getFieldData();
        assert(nn * 3 + dit->get()->getDofCoeffIdx() < 9);
        dispDofsColIdx[nn * 3 + dit->get()->getDofCoeffIdx()] =
            dit->get()->getPetscGlobalDofIdx();
      }
 
      if (trans) {
        auto dit = loIt(col_dofs, field_bit_number, conn_face[nn]);
        auto hi_dit = hiIt(col_dofs, field_bit_number, conn_face[nn]);
        if (std::distance(dit, hi_dit) > 0) {
          if (std::distance(dit, hi_dit) != 1) {
            SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                     "data inconsistency, number of dof on node for < %s > "
                     "should be 1",
                     lambdaFieldName.c_str());
          }
          if (dit->get()->getPetscLocalDofIdx() < 0) {
            SETERRQ(
                PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "data inconsistency, negative index of local dofs on element");
          }
          lambdaDofsColIdx[nn] = dit->get()->getPetscGlobalDofIdx();
        }
      }
    }
    MoFEMFunctionReturn(0);
  }
 
  Range problemTets;
 
  /**
   * \brief calculate direvatives
   *
   */
  PetscErrorCode calcDirevatives(double *diffNTRI, double *dofs_X,
                                 double *dofs_iX, double *C, double *iC,
                                 double *T, double *iT) {
    MoFEMFunctionBegin;
    double diffX_xi[3], diffX_eta[3];
    bzero(diffX_xi, 3 * sizeof(double));
    bzero(diffX_eta, 3 * sizeof(double));
    double i_diffX_xi[3], i_diffX_eta[3];
    bzero(i_diffX_xi, 3 * sizeof(double));
    bzero(i_diffX_eta, 3 * sizeof(double));
    Range adj_side_elems;
 
    BitRefLevel bit = problemPtr->getBitRefLevel();
    if (!nInTheLoop) {
      problemTets.clear();
      CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
          bit, BitRefLevel().set(), MBTET, problemTets);
    }
 
    CHKERR mField.getInterface<BitRefManager>()->getAdjacencies(
        bit, &face, 1, 3, adj_side_elems);
    adj_side_elems = intersect(problemTets, adj_side_elems);
    // adj_side_elems = adj_side_elems.subset_by_type(MBTET);
    if (adj_side_elems.size() == 0) {
      Range adj_tets_on_surface;
      BitRefLevel bit_tet_on_surface;
      bit_tet_on_surface.set(BITREFLEVEL_SIZE - 2);
      CHKERR mField.getInterface<BitRefManager>()->getAdjacencies(
          bit_tet_on_surface, &face, 1, 3, adj_tets_on_surface,
          moab::Interface::INTERSECT, 0);
 
      adj_side_elems.insert(*adj_tets_on_surface.begin());
    }
    if (adj_side_elems.size() != 1) {
      adj_side_elems.clear();
      CHKERR mField.getInterface<BitRefManager>()->getAdjacencies(
          bit, &face, 1, 3, adj_side_elems, moab::Interface::INTERSECT, 5);
 
      Range::iterator it = adj_side_elems.begin();
      for (; it != adj_side_elems.end(); it++) {
        Range nodes;
        CHKERR mField.get_moab().get_connectivity(&*it, 1, nodes, true);
        PetscPrintf(PETSC_COMM_WORLD, "%lu %lu %lu %lu\n", nodes[0], nodes[1],
                    nodes[2], nodes[3]);
      }
      ParallelComm *pcomm =
          ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
      if (pcomm->rank() == 0) {
        EntityHandle out_meshset;
        CHKERR mField.get_moab().create_meshset(MESHSET_SET, out_meshset);
        CHKERR mField.get_moab().add_entities(out_meshset, adj_side_elems);
        CHKERR mField.get_moab().add_entities(out_meshset, &face, 1);
        CHKERR mField.get_moab().write_file("debug_error.vtk", "VTK", "",
                                            &out_meshset, 1);
        CHKERR mField.get_moab().delete_entities(&out_meshset, 1);
      }
      SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "expect 1 tet but is %u", adj_side_elems.size());
    }
    EntityHandle side_elem = *adj_side_elems.begin();
    int order[] = {0, 1, 2};
    if (side_elem != 0) {
      int side_number, sense, offset;
      CHKERR moab.side_number(side_elem, face, side_number, sense, offset);
      if (sense == -1) {
        order[0] = 1;
        order[1] = 0;
      }
    }
    // calculate tangent vectors
    // those vectors are in plane of face
    for (int nn = 0; nn < 3; nn++) {
      diffX_xi[0] +=
          dofs_X[3 * order[nn] + 0] * diffNTRI[2 * nn + 0]; // unit [ m ]
      diffX_xi[1] += dofs_X[3 * order[nn] + 1] * diffNTRI[2 * nn + 0];
      diffX_xi[2] += dofs_X[3 * order[nn] + 2] * diffNTRI[2 * nn + 0];
      diffX_eta[0] += dofs_X[3 * order[nn] + 0] * diffNTRI[2 * nn + 1];
      diffX_eta[1] += dofs_X[3 * order[nn] + 1] * diffNTRI[2 * nn + 1];
      diffX_eta[2] += dofs_X[3 * order[nn] + 2] * diffNTRI[2 * nn + 1];
      if (dofs_iX == NULL)
        continue;
      i_diffX_xi[0] += dofs_iX[3 * order[nn] + 0] * diffNTRI[2 * nn + 0];
      i_diffX_xi[1] += dofs_iX[3 * order[nn] + 1] * diffNTRI[2 * nn + 0];
      i_diffX_xi[2] += dofs_iX[3 * order[nn] + 2] * diffNTRI[2 * nn + 0];
      i_diffX_eta[0] += dofs_iX[3 * order[nn] + 0] * diffNTRI[2 * nn + 1];
      i_diffX_eta[1] += dofs_iX[3 * order[nn] + 1] * diffNTRI[2 * nn + 1];
      i_diffX_eta[2] += dofs_iX[3 * order[nn] + 2] * diffNTRI[2 * nn + 1];
    }
    // spins
    double SpinX_xi[9], SpinX_eta[9];
    CHKERR Spin(SpinX_xi, diffX_xi);
    CHKERR Spin(SpinX_eta, diffX_eta);
    double iSpinX_xi[9], iSpinX_eta[9];
    CHKERR Spin(iSpinX_xi, i_diffX_xi);
    CHKERR Spin(iSpinX_eta, i_diffX_eta);
    __CLPK_doublecomplex xSpinX_xi[9], xSpinX_eta[9];
    CHKERR make_complex_matrix(SpinX_xi, iSpinX_xi, xSpinX_xi);    // unit [ m ]
    CHKERR make_complex_matrix(SpinX_eta, iSpinX_eta, xSpinX_eta); // unit [ m ]
    // calculate complex face normal vector
    __CLPK_doublecomplex x_one = {1, 0};
    __CLPK_doublecomplex x_zero = {0, 0};
    __CLPK_doublecomplex x_normal[3];
    __CLPK_doublecomplex x_diffX_eta[3];
    x_diffX_eta[0].r = diffX_eta[0];
    x_diffX_eta[0].i = i_diffX_eta[0];
    x_diffX_eta[1].r = diffX_eta[1];
    x_diffX_eta[1].i = i_diffX_eta[1];
    x_diffX_eta[2].r = diffX_eta[2];
    x_diffX_eta[2].i = i_diffX_eta[2];
    cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_one, xSpinX_xi, 3,
                x_diffX_eta, 1, &x_zero, x_normal, 1); // unit [ m^2 ]
    /*//debug
    Tag th_normal1,th_normal2;
    double def_NORMAL[] = {0, 0, 0};
    CHKERR
    moab.tag_get_handle("NORMAL1",3,MB_TYPE_DOUBLE,th_normal1,MB_TAG_CREAT|MB_TAG_SPARSE,def_NORMAL);
    CHKERR
    moab.tag_get_handle("NORMAL2",3,MB_TYPE_DOUBLE,th_normal2,MB_TAG_CREAT|MB_TAG_SPARSE,def_NORMAL);
    double normal1[3] = { x_normal[0].r, x_normal[1].r, x_normal[2].r };
    double normal2[3];
    CHKERR ShapeFaceNormalMBTRI(&diffNTRI[0],&*coords.data().begin(),normal2);
    if(side_elem!=0) {
      int side_number,sense,offset;
      CHKERR moab.side_number(side_elem,face,side_number,sense,offset);
      if(sense == -1) {
        cblas_dscal(3,-1,normal2,1);
      }
    }
    double normal1_nrm2 = cblas_dnrm2(3,normal1,1);
    cblas_dscal(3,1./normal1_nrm2,normal1,1);
    double normal2_nrm2 = cblas_dnrm2(3,normal2,1);
    cblas_dscal(3,1./normal2_nrm2,normal2,1);
    CHKERR moab.tag_set_data(th_normal1,&face,1,normal1);
    CHKERR moab.tag_set_data(th_normal2,&face,1,normal2); */
    // calculate complex normal length
    // double __complex__ x_nrm2 = csqrt(
    //   cpow((x_normal[0].r+I*x_normal[0].i),2+I*0)+
    //   cpow((x_normal[1].r+I*x_normal[1].i),2+I*0)+
    //   cpow((x_normal[2].r+I*x_normal[2].i),2+I*0)
    // );
    complex<double> x_nrm2;
    x_nrm2 = sqrt(pow(complex<double>(x_normal[0].r, x_normal[0].i), 2) +
                  pow(complex<double>(x_normal[1].r, x_normal[1].i), 2) +
                  pow(complex<double>(x_normal[2].r, x_normal[2].i), 2));
    // calculate dA/dX
    __CLPK_doublecomplex xNSpinX_xi[3], xNSpinX_eta[3];
    bzero(xNSpinX_xi, 3 * sizeof(__CLPK_doublecomplex));
    bzero(xNSpinX_eta, 3 * sizeof(__CLPK_doublecomplex));
    // __CLPK_doublecomplex x_scalar = { -creal(1./x_nrm2), -cimag(1./x_nrm2) };
    // // unit [ 1/m^2 ]
    __CLPK_doublecomplex x_scalar = {-real(1. / x_nrm2),
                                     -imag(1. / x_nrm2)}; // unit [ 1/m^2 ]
    cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_scalar, xSpinX_xi, 3,
                x_normal, 1, &x_zero, xNSpinX_xi,
                1); // unit [ 1/m^2 * m = 1/m ]
    cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_scalar, xSpinX_eta, 3,
                x_normal, 1, &x_zero, xNSpinX_eta,
                1); // unit [ 1/m^2 * m = 1/m ]
    if (C != NULL)
      bzero(C, 9 * sizeof(double));
    if (iC != NULL)
      bzero(iC, 9 * sizeof(double));
    if (T != NULL)
      bzero(T, 9 * sizeof(double));
    if (iT != NULL)
      bzero(iT, 9 * sizeof(double));
    for (int nn = 0; nn < 3; nn++) {
      double A[3], iA[3];
      A[0] = xNSpinX_xi[0].r * diffNTRI[2 * order[nn] + 1] -
             xNSpinX_eta[0].r * diffNTRI[2 * order[nn] + 0]; // unit [ 1/m ]
      A[1] = xNSpinX_xi[1].r * diffNTRI[2 * order[nn] + 1] -
             xNSpinX_eta[1].r * diffNTRI[2 * order[nn] + 0];
      A[2] = xNSpinX_xi[2].r * diffNTRI[2 * order[nn] + 1] -
             xNSpinX_eta[2].r * diffNTRI[2 * order[nn] + 0];
      iA[0] = xNSpinX_xi[0].i * diffNTRI[2 * order[nn] + 1] -
              xNSpinX_eta[0].i * diffNTRI[2 * order[nn] + 0]; // unit [ 1/m ]
      iA[1] = xNSpinX_xi[1].i * diffNTRI[2 * order[nn] + 1] -
              xNSpinX_eta[1].i * diffNTRI[2 * order[nn] + 0];
      iA[2] = xNSpinX_xi[2].i * diffNTRI[2 * order[nn] + 1] -
              xNSpinX_eta[2].i * diffNTRI[2 * order[nn] + 0];
      if (C != NULL) {
        C[3 * nn + 0] = A[0];
        C[3 * nn + 1] = A[1];
        C[3 * nn + 2] = A[2];
      }
      if (iC != NULL) {
        iC[3 * nn + 0] = iA[0];
        iC[3 * nn + 1] = iA[1];
        iC[3 * nn + 2] = iA[2];
      }
      if ((T != NULL) || (iT != NULL)) {
        double SpinA[9];
        CHKERR Spin(SpinA, A); // unit [1/m]
        double iSpinA[9];
        CHKERR Spin(iSpinA, iA); // unit [1/m]
        __CLPK_doublecomplex xSpinA[9];
        // make spin matrix to calculate cross product
        CHKERR make_complex_matrix(SpinA, iSpinA, xSpinA);
        __CLPK_doublecomplex xT[3];
        cblas_zgemv(CblasRowMajor, CblasNoTrans, 3, 3, &x_scalar, xSpinA, 3,
                    x_normal, 1, &x_zero, xT, 1); // unit [1/m]
        if (T != NULL) {
          T[3 * nn + 0] = xT[0].r;
          T[3 * nn + 1] = xT[1].r;
          T[3 * nn + 2] = xT[2].r;
        }
        if (iT != NULL) {
          iT[3 * nn + 0] = xT[0].i;
          iT[3 * nn + 1] = xT[1].i;
          iT[3 * nn + 2] = xT[2].i;
        }
      }
    }
    if (C != NULL)
      cblas_dscal(9, 0.25, C, 1);
    if (iC != NULL)
      cblas_dscal(9, 0.25, iC, 1);
    if (T != NULL)
      cblas_dscal(9, 0.25, T, 1);
    if (iT != NULL)
      cblas_dscal(9, 0.25, iT, 1);
    MoFEMFunctionReturn(0);
  }
 
  map<DofIdx, Vec> mapV;
  PetscErrorCode preProcess() {
    MoFEMFunctionBegin;
 
    if (Q != PETSC_NULL) {
      const auto field_bit_number =
          mField.get_field_bit_number(lambdaFieldName);
      for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(
               problemPtr, field_bit_number, dofs_it)) {
        CHKERR MatCreateVecs(C, &mapV[dofs_it->get()->getPetscGlobalDofIdx()],
                             PETSC_NULL);
        CHKERR VecZeroEntries(mapV[dofs_it->get()->getPetscGlobalDofIdx()]);
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
  PetscErrorCode operator()() {
    MoFEMFunctionBegin;
 
    CHKERR getData(true, false);
 
    ublas::vector<double, ublas::bounded_array<double, 9>> ELEM_CONSTRAIN(9);
    CHKERR calcDirevatives(&*diffNTRI.data().begin(), &*dofs_X.data().begin(),
                           NULL, &*ELEM_CONSTRAIN.data().begin(), NULL, NULL,
                           NULL);
    for (int nn = 0; nn < 3; nn++) {
      if (lambdaDofsRowIdx[nn] == -1)
        continue;
      for (int NN = 0; NN < 3; NN++) {
        if (NN != nn)
          continue;
        if (Q == PETSC_NULL) {
          CHKERR MatSetValues(C, 1, &(lambdaDofsRowIdx.data()[nn]), 3,
                              &(dispDofsColIdx.data()[3 * NN]),
                              &ELEM_CONSTRAIN.data()[3 * NN], ADD_VALUES);
        }
        if (Q != PETSC_NULL) {
          if (mapV.find(lambdaDofsRowIdx[nn]) == mapV.end())
            SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
          CHKERR VecSetValues(mapV[lambdaDofsRowIdx[nn]], 3,
                              &(dispDofsColIdx.data()[3 * NN]),
                              &ELEM_CONSTRAIN.data()[3 * NN], ADD_VALUES);
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
  PetscErrorCode postProcess() {
    MoFEMFunctionBegin;
    if (Q != PETSC_NULL) {
      ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
      Vec Qv;
      CHKERR mField.getInterface<VecManager>()->vecCreateGhost(
          problemPtr->getName(), COL, &Qv);
      const auto field_bit_number =
          mField.get_field_bit_number(lambdaFieldName);
      for (_IT_NUMEREDDOF_ROW_BY_BITNUMBER_FOR_LOOP_(problemPtr,
                                                     field_bit_number, dofs)) {
        if (mapV.find(dofs->get()->getPetscGlobalDofIdx()) == mapV.end()) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "data inconsistency");
        }
        CHKERR VecAssemblyBegin(mapV[dofs->get()->getPetscGlobalDofIdx()]);
        CHKERR VecAssemblyEnd(mapV[dofs->get()->getPetscGlobalDofIdx()]);
        CHKERR MatMult(Q, mapV[dofs->get()->getPetscGlobalDofIdx()], Qv);
        CHKERR VecGhostUpdateBegin(Qv, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(Qv, INSERT_VALUES, SCATTER_FORWARD);
        double *array;
        CHKERR VecGetArray(Qv, &array);
        vector<DofIdx> glob_idx;
        vector<double> vals;
        for (_IT_NUMEREDDOF_COL_BY_ENT_FOR_LOOP_(problemPtr,
                                                 dofs->get()->getEnt(), diit)) {
          if (diit->get()->getPart() != pcomm->rank())
            continue;
          if (diit->get()->getName() != "MESH_NODE_POSITIONS")
            continue;
          glob_idx.push_back(diit->get()->getPetscGlobalDofIdx());
          if (diit->get()->getPetscLocalDofIdx() < 0) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "data inconsistency");
          }
          vals.push_back(array[diit->get()->getPetscLocalDofIdx()]);
        }
        int row = dofs->get()->getPetscGlobalDofIdx();
        CHKERR MatSetValues(C, 1, &row, glob_idx.size(),
                            glob_idx.empty() ? PETSC_NULL : &glob_idx[0],
                            vals.empty() ? PETSC_NULL : &vals[0], ADD_VALUES);
        CHKERR VecRestoreArray(Qv, &array);
        CHKERR VecDestroy(&mapV[dofs->get()->getPetscGlobalDofIdx()]);
      }
      CHKERR VecDestroy(&Qv);
      mapV.clear();
    }
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * Constrains is vector_tangent_to_front*nodal_element_quality_at_crack_front =
 * 0, multiplying that constrain by Legrange multipliers, we get
 *
 * 1) vector_tangent_to_front*nodal_element_quality_at_crack_front=0 for each
 * lambda, 2) lambda*vector_tangent_to_front*direvatives_of_quality = 0
 *
 * Those two equations need to be linearsed and assembled into system of
 * equations
 *
 *
 */
struct TangentWithMeshSmoothingFrontConstrain : public ConstrainConstantAarea {
 
  Vec frontF;
  Vec tangentFrontF;
  const double eps;
  bool ownVectors;
 
  TangentWithMeshSmoothingFrontConstrain(MoFEM::Interface &m_field,
                                         string lambda_field_name,
                                         int verbose = 0)
      : ConstrainConstantAarea(m_field, PETSC_NULL, PETSC_NULL,
                               lambda_field_name, verbose),
        frontF(PETSC_NULL), tangentFrontF(PETSC_NULL), eps(1e-10),
        ownVectors(false) {}
 
  virtual ~TangentWithMeshSmoothingFrontConstrain() {
    if (ownVectors) {
      CHKERR VecDestroy(&frontF);
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
      CHKERR VecDestroy(&tangentFrontF);
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
  }
 
  Tag thFrontTangent;
  PetscErrorCode preProcess() {
    MoFEMFunctionBegin;
 
    switch (ts_ctx) {
    case CTX_TSSETIFUNCTION: {
      snes_ctx = CTX_SNESSETFUNCTION;
      snes_f = ts_F;
      break;
    }
    case CTX_TSSETIJACOBIAN: {
      snes_ctx = CTX_SNESSETJACOBIAN;
      snes_B = ts_B;
      break;
    }
    default:
      break;
    }
 
    CHKERR ConstrainConstantAarea::preProcess();
 
    switch (snes_ctx) {
    case CTX_SNESSETFUNCTION: {
      CHKERR VecAssemblyBegin(snes_f);
      CHKERR VecAssemblyEnd(snes_f);
      CHKERR VecZeroEntries(tangentFrontF);
      CHKERR VecGhostUpdateBegin(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecAssemblyBegin(tangentFrontF);
      CHKERR VecAssemblyEnd(tangentFrontF);
    } break;
    case CTX_SNESSETJACOBIAN: {
      CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
      CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
    } break;
    default: {
    } break;
    }
    MoFEMFunctionReturn(0);
  }
 
  PetscErrorCode operator()() {
    MoFEMFunctionBegin;
    // get elemet hanlde
    EntityHandle face = numeredEntFiniteElementPtr->getEnt();
    // get adjacent tets
    Range tet;
    CHKERR mField.get_moab().get_adjacencies(&face, 1, 3, false, tet);
    // get data and indices on dofs
    CHKERR getData(true, true);
 
    // Calculate tangent vector constrains
    ublas::vector<double, ublas::bounded_array<double, 9>> tangent_constrains(
        9);
    CHKERR calcDirevatives(&*diffNTRI.data().begin(), &*dofs_X.data().begin(),
                           NULL, NULL, NULL,
                           &*tangent_constrains.data().begin(), NULL);
    // take in account face orientation in respect crack surface
    Tag th_interface_side;
    CHKERR moab.tag_get_handle("INTERFACE_SIDE", th_interface_side);
    int side;
    CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
    if (side == 1) {
      tangent_constrains *= -1;
    }
    // get smoothing element forces
    ublas::vector<double, ublas::bounded_array<double, 9>>
        f_front_mesh_smoothing(9);
    ublas::noalias(f_front_mesh_smoothing) = ublas::zero_vector<double>(9);
    double *f_front_mesh_array;
    CHKERR VecGetArray(frontF, &f_front_mesh_array);
    for (int nn = 0; nn < 3; nn++) {
      for (int dd = 0; dd < 3; dd++) {
        if (localDispDofsRowIdx[3 * nn + dd] == -1)
          continue;
        f_front_mesh_smoothing[3 * nn + dd] =
            f_front_mesh_array[localDispDofsRowIdx[3 * nn + dd]];
      }
    }
    CHKERR VecRestoreArray(frontF, &f_front_mesh_array);
    // calculate tangent matrix
    if (snes_ctx == CTX_SNESSETJACOBIAN) {
      // get radius
      double center[3];
      tricircumcenter3d_tp(&coords.data()[0], &coords.data()[3],
                           &coords.data()[6], center, NULL, NULL);
      cblas_daxpy(3, -1, &coords.data()[0], 1, center, 1);
      double r = cblas_dnrm2(3, center, 1);
      for (int nn = 0; nn < 3; nn++) {
        for (int dd = 0; dd < 3; dd++) {
          // ---> calculate tangent starts here
          ublas::vector<double, ublas::bounded_array<double, 9>> idofs_X(9, 0);
          idofs_X[nn * 3 + dd] = r * eps;
          ublas::vector<double, ublas::bounded_array<double, 9>>
              direvative_of_constrains(9);
          // use complex derivative
          CHKERR calcDirevatives(&*diffNTRI.data().begin(),
                                 &*dofs_X.data().begin(),
                                 &*idofs_X.data().begin(), NULL, NULL, NULL,
                                 &*direvative_of_constrains.data().begin());
          if (side == 1) {
            direvative_of_constrains /= -r * eps;
          } else {
            direvative_of_constrains /= +r * eps;
          }
          // dg -> C*q_quality
          double g[3] = {0, 0, 0};
          for (int nnn = 0; nnn < 3; nnn++) {
            if (lambdaDofsRowIdx[nnn] == -1)
              continue;
            for (int ddd = 0; ddd < 3; ddd++) {
              g[nnn] += direvative_of_constrains[nnn * 3 + ddd] *
                        f_front_mesh_smoothing[3 * nnn + ddd];
            }
          }
          for (int nnn = 0; nnn < 3; nnn++) {
            if (lambdaDofsRowIdx[nnn] == -1)
              continue;
            CHKERR MatSetValues(snes_B, 1, &lambdaDofsRowIdx[nnn], 1,
                                &dispDofsColIdx[3 * nn + dd], &g[nnn],
                                ADD_VALUES);
          }
          // dCT_lambda
          for (int nnn = 0; nnn < 3; nnn++) {
            for (int ddd = 0; ddd < 3; ddd++) {
              direvative_of_constrains[nnn * 3 + ddd] *= lambda[nnn];
            }
          }
          for (int nnn = 0; nnn < 3; nnn++) {
            if (lambdaDofsRowIdx[nnn] == -1)
              continue;
            CHKERR MatSetValues(snes_B, 3, &dispDofsRowIdx[3 * nnn], 1,
                                &dispDofsColIdx[3 * nn + dd],
                                &direvative_of_constrains[3 * nnn], ADD_VALUES);
          }
          // ---> calculate tangent end here
        }
      }
    }
    switch (snes_ctx) {
    case CTX_SNESSETFUNCTION: {
      ublas::vector<double, ublas::bounded_array<double, 3>> g(3);
      for (int nn = 0; nn < 3; nn++) {
        g[nn] = cblas_ddot(3, &tangent_constrains[3 * nn], 1,
                           &f_front_mesh_smoothing[3 * nn], 1);
      }
      // cerr << "g : " << g << endl;
      CHKERR VecSetValues(snes_f, 3, &*lambdaDofsRowIdx.data().begin(),
                          &*g.data().begin(), ADD_VALUES);
      CHKERR VecSetValues(tangentFrontF, 9, &dispDofsRowIdx[0],
                          &*tangent_constrains.data().begin(), ADD_VALUES);
      ublas::vector<double, ublas::bounded_array<double, 9>> f(9);
      for (int nn = 0; nn < 3; nn++) {
        for (int dd = 0; dd < 3; dd++) {
          f[nn * 3 + dd] = lambda[nn] * tangent_constrains[3 * nn + dd];
        }
      }
      // cerr << "f : " << f << endl;
      CHKERR VecSetValues(snes_f, 9, &dispDofsRowIdx[0], &*f.data().begin(),
                          ADD_VALUES);
      /*//TAG - only for one proc analysis
      for(int nn = 0;nn<3;nn++) {
      EntityHandle ent = lambda_dofs_row_ents[nn];
      if(ent == no_handle) continue;
      double *t;
      CHKERR mField.get_moab().tag_get_by_ptr(thFrontTangent,&ent,1,(const void
    **)&t); cblas_daxpy(3,+1,&tangent_constrains[3*nn],1,t,1);
    }*/
    } break;
    case CTX_SNESSETJACOBIAN: {
      for (int nn = 0; nn < 3; nn++) {
        int lambda_dof_idx = lambdaDofsColIdx[nn];
        CHKERR MatSetValues(snes_B, 3, &dispDofsRowIdx[3 * nn], 1,
                            &lambda_dof_idx, &tangent_constrains[3 * nn],
                            ADD_VALUES);
      }
    } break;
    default:
      break;
    }
    MoFEMFunctionReturn(0);
  }
 
  PetscErrorCode postProcess() {
    MoFEMFunctionBegin;
    switch (snes_ctx) {
    case CTX_SNESSETFUNCTION: {
      CHKERR VecAssemblyBegin(tangentFrontF);
      CHKERR VecAssemblyEnd(tangentFrontF);
      CHKERR VecGhostUpdateBegin(tangentFrontF, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(tangentFrontF, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateBegin(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(tangentFrontF, INSERT_VALUES, SCATTER_FORWARD);
    } break;
    case CTX_SNESSETJACOBIAN: {
      CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
      CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
    } break;
    default: {
    } break;
    }
    MoFEMFunctionReturn(0);
  }
};
 
} // namespace ObosleteUsersModules
 
#endif // __COMPLEX_CONST_AREA_HPP__