ObosleteUsersModules::ConstrainConstantAarea Struct Reference

#include <users_modules/fracture_mechanics/src/ComplexConstArea.hpp>

Inheritance diagram for ObosleteUsersModules::ConstrainConstantAarea:

Collaboration diagram for ObosleteUsersModules::ConstrainConstantAarea:

Public Member Functions
	ConstrainConstantAarea (MoFEM::Interface &m_field, Mat _C, Mat _Q, string lambda_field_name, int _verbose=0)
virtual	~ConstrainConstantAarea ()=default
template<typename DOFS >
auto	loIt (DOFS dofs, const FieldBitNumber bit, const EntityHandle ent)
template<typename DOFS >
auto	hiIt (DOFS dofs, const FieldBitNumber bit, const EntityHandle ent)
PetscErrorCode	getData (bool is_that_C_otherwise_dC, bool trans)
	get face data, indices and coors and nodal values More...
PetscErrorCode	calcDirevatives (double *diffNTRI, double *dofs_X, double *dofs_iX, double *C, double *iC, double *T, double *iT)
	calculate direvatives More...
PetscErrorCode	preProcess ()
PetscErrorCode	operator() ()
PetscErrorCode	postProcess ()

Public Attributes
MoFEM::Interface &	mField
moab::Interface &	moab
Mat	C
Mat	Q
string	lambdaFieldName
int	verbose
MatrixDouble	diffNTRI
const double *	G_TRI_W
vector< DofIdx >	DirichletBC
VectorInt	dispDofsColIdx
VectorInt	dispDofsRowIdx
VectorInt	localDispDofsRowIdx
VectorInt	lambdaDofsRowIdx
VectorInt	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
EntityHandle	face
Range	problemTets
map< DofIdx, Vec >	mapV

Detailed Description

dN/dX = (1/A) * [ Spin[dX/dksi]*dN/deta - Spin[dX/deta]*dN/dksi ]

Definition at line 37 of file ComplexConstArea.hpp.

Constructor & Destructor Documentation

ConstrainConstantAarea()

ObosleteUsersModules::ConstrainConstantAarea::ConstrainConstantAarea ( MoFEM::Interface &  m_field, Mat  _C, Mat  _Q, string  lambda_field_name, int  _verbose = 0  )

inline

Definition at line 48 of file ComplexConstArea.hpp.

      : 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);
  }

~ConstrainConstantAarea()

virtual ObosleteUsersModules::ConstrainConstantAarea::~ConstrainConstantAarea ( )

virtualdefault

Member Function Documentation

calcDirevatives()

PetscErrorCode ObosleteUsersModules::ConstrainConstantAarea::calcDirevatives ( double *  diffNTRI, double *  dofs_X, double *  dofs_iX, double *  C, double *  iC, double *  T, double *  iT  )

inline

calculate direvatives

Definition at line 245 of file ComplexConstArea.hpp.

                                                        {
    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);
  }

getData()

PetscErrorCode ObosleteUsersModules::ConstrainConstantAarea::getData ( bool  is_that_C_otherwise_dC, bool  trans  )

inline

get face data, indices and coors and nodal values

Parameters

is_that_C_otherwise_dC

Definition at line 105 of file ComplexConstArea.hpp.

                                                                  {
    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);
  }

hiIt()

template<typename DOFS >

auto ObosleteUsersModules::ConstrainConstantAarea::hiIt ( DOFS  dofs, const FieldBitNumber  bit, const EntityHandle  ent  )

inline

Definition at line 95 of file ComplexConstArea.hpp.

                                           {
    return dofs->upper_bound(DofEntity::getHiFieldEntityUId(bit, ent));
  };

loIt()

template<typename DOFS >

auto ObosleteUsersModules::ConstrainConstantAarea::loIt ( DOFS  dofs, const FieldBitNumber  bit, const EntityHandle  ent  )

inline

Definition at line 89 of file ComplexConstArea.hpp.

                                           {
    return dofs->lower_bound(DofEntity::getLoFieldEntityUId(bit, ent));
  };

operator()()

PetscErrorCode ObosleteUsersModules::ConstrainConstantAarea::operator() ( )

inline

Definition at line 487 of file ComplexConstArea.hpp.

                              {
    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);
  }

postProcess()

PetscErrorCode ObosleteUsersModules::ConstrainConstantAarea::postProcess ( )

inline

Definition at line 520 of file ComplexConstArea.hpp.

                               {
    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);
  }

preProcess()

PetscErrorCode ObosleteUsersModules::ConstrainConstantAarea::preProcess ( )

inline

Definition at line 470 of file ComplexConstArea.hpp.

                              {
    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);
  }

Member Data Documentation

C

Mat ObosleteUsersModules::ConstrainConstantAarea::C

Definition at line 43 of file ComplexConstArea.hpp.

coords

ublas::vector<double, ublas::bounded_array<double, 9> > ObosleteUsersModules::ConstrainConstantAarea::coords

Definition at line 82 of file ComplexConstArea.hpp.

diffNTRI

MatrixDouble ObosleteUsersModules::ConstrainConstantAarea::diffNTRI

Definition at line 76 of file ComplexConstArea.hpp.

DirichletBC

vector<DofIdx> ObosleteUsersModules::ConstrainConstantAarea::DirichletBC

Definition at line 78 of file ComplexConstArea.hpp.

dispDofsColIdx

VectorInt ObosleteUsersModules::ConstrainConstantAarea::dispDofsColIdx

Definition at line 79 of file ComplexConstArea.hpp.

dispDofsRowIdx

VectorInt ObosleteUsersModules::ConstrainConstantAarea::dispDofsRowIdx

Definition at line 79 of file ComplexConstArea.hpp.

dofs_X

ublas::vector<double, ublas::bounded_array<double, 9> > ObosleteUsersModules::ConstrainConstantAarea::dofs_X

Definition at line 83 of file ComplexConstArea.hpp.

face

EntityHandle ObosleteUsersModules::ConstrainConstantAarea::face

Definition at line 86 of file ComplexConstArea.hpp.

G_TRI_W

const double* ObosleteUsersModules::ConstrainConstantAarea::G_TRI_W

Definition at line 77 of file ComplexConstArea.hpp.

lambda

ublas::vector<double, ublas::bounded_array<double, 3> > ObosleteUsersModules::ConstrainConstantAarea::lambda

Definition at line 84 of file ComplexConstArea.hpp.

lambdaDofsColIdx

VectorInt ObosleteUsersModules::ConstrainConstantAarea::lambdaDofsColIdx

Definition at line 81 of file ComplexConstArea.hpp.

lambdaDofsRowIdx

VectorInt ObosleteUsersModules::ConstrainConstantAarea::lambdaDofsRowIdx

Definition at line 81 of file ComplexConstArea.hpp.

lambdaFieldName

string ObosleteUsersModules::ConstrainConstantAarea::lambdaFieldName

Definition at line 45 of file ComplexConstArea.hpp.

localDispDofsRowIdx

VectorInt ObosleteUsersModules::ConstrainConstantAarea::localDispDofsRowIdx

Definition at line 80 of file ComplexConstArea.hpp.

mapV

map<DofIdx, Vec> ObosleteUsersModules::ConstrainConstantAarea::mapV

Definition at line 469 of file ComplexConstArea.hpp.

mField

MoFEM::Interface& ObosleteUsersModules::ConstrainConstantAarea::mField

Definition at line 39 of file ComplexConstArea.hpp.

moab

moab::Interface& ObosleteUsersModules::ConstrainConstantAarea::moab

Definition at line 41 of file ComplexConstArea.hpp.

problemTets

Range ObosleteUsersModules::ConstrainConstantAarea::problemTets

Definition at line 239 of file ComplexConstArea.hpp.

Q

Mat ObosleteUsersModules::ConstrainConstantAarea::Q

Definition at line 44 of file ComplexConstArea.hpp.

verbose

int ObosleteUsersModules::ConstrainConstantAarea::verbose

Definition at line 46 of file ComplexConstArea.hpp.

---

The documentation for this struct was generated from the following file:

• users_modules/fracture_mechanics/src/ComplexConstArea.hpp