Operator the left hand side matrix. More...

Inheritance diagram for OpLhsSkeleton:

Collaboration diagram for OpLhsSkeleton:

Public Member Functions
	OpLhsSkeleton (boost::shared_ptr< FaceSideEle > side_fe_ptr, boost::shared_ptr< MatrixDouble > d_mat_fluid_ptr, boost::shared_ptr< MatrixDouble > d_mat_solid_ptr, boost::shared_ptr< MatrixDouble > d_mat_visc_fluid_ptr, boost::shared_ptr< MatrixDouble > d_mat_visc_solid_ptr, boost::shared_ptr< Range > r_visc_fluid_ptr, boost::shared_ptr< Range > r_visc_solid_ptr, SmartPetscObj< Mat > K, SmartPetscObj< Mat > C, SmartPetscObj< Mat > M)
	Construct a new OpH1LhsSkeleton. More...

MoFEMErrorCode	doWork (int side, EntityType type, EntitiesFieldData::EntData &data)

Private Attributes
boost::shared_ptr< FaceSideEle >	sideFEPtr
	pointer to element to get data on edge/face sides More...

MatrixDouble	locMatK

MatrixDouble	locMatC

MatrixDouble	locMatM

boost::shared_ptr< MatrixDouble >	dMatFluidPtr

boost::shared_ptr< MatrixDouble >	dMatSolidPtr

boost::shared_ptr< MatrixDouble >	dViscFluidPtr

boost::shared_ptr< MatrixDouble >	dViscSolidPtr

boost::shared_ptr< Range >	rFluidPtr

boost::shared_ptr< Range >	rSolidPtr

SmartPetscObj< Mat >	K

SmartPetscObj< Mat >	C

SmartPetscObj< Mat >	M

Detailed Description

Operator the left hand side matrix.

Examples: level_set.cpp.

Definition at line 110 of file fluid_structure_eigenproblem.cpp.

Constructor & Destructor Documentation

◆ OpLhsSkeleton()

OpLhsSkeleton::OpLhsSkeleton	(	boost::shared_ptr< FaceSideEle >	side_fe_ptr,
		boost::shared_ptr< MatrixDouble >	mat_d_fluid_ptr,
		boost::shared_ptr< MatrixDouble >	mat_d_solid_ptr,
		boost::shared_ptr< MatrixDouble >	mat_d_visc_fluid_ptr,
		boost::shared_ptr< MatrixDouble >	mat_d_visc_solid_ptr,
		boost::shared_ptr< Range >	r_fluid_ptr,
		boost::shared_ptr< Range >	r_solid_ptr,
		SmartPetscObj< Mat >	K,
		SmartPetscObj< Mat >	C,
		SmartPetscObj< Mat >	M
	)

Construct a new OpH1LhsSkeleton.

Parameters

side_fe_ptr pointer to FE to evaluate side elements

Parameters

side_fe_ptr pointer to FE to evaluate side elements

Definition at line 1042 of file fluid_structure_eigenproblem.cpp.

    : SkeletonEleOp(NOSPACE, SkeletonEleOp::OPSPACE), sideFEPtr(side_fe_ptr),
      dMatFluidPtr(mat_d_fluid_ptr), dMatSolidPtr(mat_d_solid_ptr),
      dViscFluidPtr(mat_d_visc_fluid_ptr), dViscSolidPtr(mat_d_visc_solid_ptr),
      rFluidPtr(r_fluid_ptr), rSolidPtr(r_solid_ptr), K(K), C(C), M(M) {}

Member Function Documentation

◆ doWork()

MoFEMErrorCode OpLhsSkeleton::doWork	(	int	side,
		EntityType	type,
		EntitiesFieldData::EntData &	data
	)

Definition at line 1055 of file fluid_structure_eigenproblem.cpp.

                                                                     {
  MoFEMFunctionBegin;
 
  // get normal of the face or edge
  auto t_normal = getFTensor1Normal();
  t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
 
  // Collect data from side domian elements
  CHKERR loopSide("dFE", sideFEPtr.get(), SPACE_DIM);
  const auto in_the_loop =
      sideFEPtr->nInTheLoop; // return number of elements on the side
 
  std::array<boost::shared_ptr<MatrixDouble>, 2> d_ptr;
  std::array<boost::shared_ptr<MatrixDouble>, 2> d_vis_ptr;
 
  std::array<std::vector<VectorInt>, 2> indices_side_map;
  std::array<std::vector<MatrixDouble>, 2> base_side_map;
  std::array<std::vector<MatrixDouble>, 2> diff_base_side_map;
 
  // MOFEM_LOG("SELF", Sev::inform)
  //     << "AAA " << sideHandle[0] << " " << sideHandle[1];
 
  // auto save_range = [&](const std::string name, const Range r) {
  //   MoFEMFunctionBegin;
  //   EntityHandle out_meshset;
  //   CHKERR getPtrFE()->mField.get_moab().create_meshset(MESHSET_SET,
  //                                                       out_meshset);
  //   CHKERR getPtrFE()->mField.get_moab().add_entities(out_meshset, r);
  //   CHKERR getPtrFE()->mField.get_moab().write_file(name.c_str(), "VTK", "",
  //                                                   &out_meshset, 1);
  //   CHKERR getPtrFE()->mField.get_moab().delete_entities(&out_meshset, 1);
  //   MoFEMFunctionReturn(0);
  // };
 
  // Tag th_n;
  // double def[] = {0, 0, 0};
  // CHKERR getPtrFE()->mField.get_moab().tag_get_handle(
  //     "NORMAL", 3, MB_TYPE_DOUBLE, th_n, MB_TAG_CREAT | MB_TAG_SPARSE, def);
 
  auto set_row_col = [&](auto a, auto b) {
    d_ptr[a] = dMatSolidPtr;
    d_ptr[b] = dMatFluidPtr;
    d_vis_ptr[a] = dViscSolidPtr;
    d_vis_ptr[b] = dViscFluidPtr;
    indices_side_map[a].swap(indicesUsSideMap[a]);
    indices_side_map[b].swap(indicesUfSideMap[b]);
    base_side_map[a].swap(baseUsSideMap[a]);
    base_side_map[b].swap(baseUfSideMap[b]);
    diff_base_side_map[a].swap(diffUsBaseSideMap[a]);
    diff_base_side_map[b].swap(diffUfBaseSideMap[b]);
  };
 
  if (rSolidPtr->find(sideHandle[0]) != rSolidPtr->end()) {
 
    if (rFluidPtr->find(sideHandle[1]) == rFluidPtr->end()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Other side should be fluid");
    }
 
    set_row_col(0, 1);
 
    if (senseMap[0] < 0)
      t_normal(i) *= -1;
 
  } else if (rFluidPtr->find(sideHandle[0]) != rFluidPtr->end()) {
 
    if (rSolidPtr->find(sideHandle[1]) == rSolidPtr->end()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Other side should be solid");
    }
 
    set_row_col(1, 0);
 
    if (senseMap[1] < 0)
      t_normal(i) *= -1;
 
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Side is not fluid or solid");
  }
 
  // auto fe_h = getFEEntityHandle();
  // CHKERR getPtrFE()->mField.get_moab().tag_set_data(th_n, &fe_h, 1,
  //                                                   &t_normal(0));
 
  // calculate  penalty
  const double s = getMeasure() / (areaMap[0] + areaMap[1]);
  const double p = pow(penalty, order) * s;
 
  // get number of integration points on face
  const size_t nb_integration_pts = getGaussPts().size2();
 
  // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
  const double beta = 1. / (in_the_loop + 1);
 
  auto integrate = [&](auto sense_row, auto &row_ind, auto &row, auto &row_diff,
                       auto sense_col, auto &col_ind, auto &col, auto &col_diff,
                       auto &d_row, auto &d_col, auto &d_vis_row,
                       auto &d_vis_col) {
    MoFEMFunctionBeginHot;
 
    // number of dofs, for homogenous approximation this value is
    // constant.
    const auto nb_rows = row_ind.size();
    const auto nb_cols = col_ind.size();
 
    const auto nb_row_base_functions = row.size2();
 
    auto t_D_row = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_row);
    auto t_D_col = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_col);
 
    auto t_vis_D_row =
        getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_row);
    auto t_vis_D_col =
        getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_col);
 
    if (nb_cols) {
 
      // resize local element matrix
      locMatK.resize(nb_rows, nb_cols, false);
      locMatK.clear();
      locMatC.resize(nb_rows, nb_cols, false);
      locMatC.clear();
      locMatM.resize(nb_rows, nb_cols, false);
      locMatM.clear();
 
      // get base functions, and integration weights
      auto t_row_base = get_ntensor(row);
      auto t_diff_row_base = get_diff_ntensor(row_diff);
 
      constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_P;
      t_P(i, j) = t_normal(i) * t_normal(j);
 
      auto t_w = getFTensor0IntegrationWeight();
 
      // iterate integration points on face/edge
      for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
        // t_w is integration weight, and measure is element area, or
        // volume, depending if problem is in 2d/3d.
        const double alpha = getMeasure() * t_w;
 
        // iterate rows
        size_t rr = 0;
        for (; rr != nb_rows / SPACE_DIM; ++rr) {
 
          auto get_vn_plus = [&](auto &t_D) {
            FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_vn_plus;
            t_vn_plus(i, j, k) =
                (beta * phi / p) *
                (t_kd(i, m) * t_D(m, j, k, l) * t_diff_row_base(l));
            return t_vn_plus;
          };
 
          auto get_vn = [&]() {
            FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_vn;
            t_vn(i, j, k) = t_normal(j) * (t_kd(i, k) * t_row_base * sense_row);
            return t_vn;
          };
 
          auto t_vn_plus = get_vn_plus(t_D_row);
          auto t_vn = get_vn();
          auto t_vn_lambda_plus = get_vn_plus(t_vis_D_row);
 
          // get base functions on columns
          auto t_col_base = get_ntensor(col, gg, 0);
          auto t_diff_col_base = get_diff_ntensor(col_diff, gg, 0);
 
          auto t_mat_K = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMatK, SPACE_DIM * rr);
          auto t_mat_C = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMatC, SPACE_DIM * rr);
          auto t_mat_M = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMatM, SPACE_DIM * rr);
 
          // iterate columns
          for (size_t cc = 0; cc != nb_cols / SPACE_DIM; ++cc) {
 
            auto get_un_plus = [&](auto &t_D) {
              FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM>
                  t_un_plus;
              t_un_plus(i, j, k) =
                  beta * ((t_P(i, m) * t_D(m, j, k, l)) * t_diff_col_base(l));
              return t_un_plus;
            };
 
            auto get_un = [&]() {
              FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_un;
              t_un(i, j, k) =
                  (t_normal(j) * (t_P(i, k) * t_col_base * sense_col));
              return t_un;
            };
 
            auto t_un_plus = get_un_plus(t_D_col);
            auto t_un = get_un();
            auto t_un_lambda_plus = get_un_plus(t_vis_D_col);
 
            // K
            t_mat_K(i, j) -=
                alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
                ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
            t_mat_K(i, j) -= alpha * (t_vn_plus(m, n, i) * t_un_plus(m, n, j));
 
            // C1
            t_mat_C(i, j) -=
                alpha * (-t_vn_lambda_plus(m, n, i)) *
                ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
            t_mat_C(i, j) -=
                alpha * (t_vn_lambda_plus(m, n, i) * t_un_plus(m, n, j));
 
            // C2
            t_mat_C(i, j) -= alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
                             ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
            t_mat_C(i, j) -=
                alpha * (t_vn_plus(m, n, i) * t_un_lambda_plus(m, n, j));
 
            // M
            t_mat_M(i, j) -= alpha * (-t_vn_lambda_plus(m, n, i)) *
                             ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
            t_mat_M(i, j) -=
                alpha * (t_vn_lambda_plus(m, n, i) * t_un_lambda_plus(m, n, j));
 
            // move to next column base and element of matrix
            ++t_col_base;
            ++t_diff_col_base;
            ++t_mat_K;
            ++t_mat_C;
            ++t_mat_M;
          }
 
          // move to next row base
          ++t_row_base;
          ++t_diff_row_base;
        }
 
        for (; rr < nb_row_base_functions; ++rr) {
          ++t_row_base;
          ++t_diff_row_base;
        }
 
        ++t_w;
      }
 
      // assemble system
 
      CHKERR ::MatSetValues(K, nb_rows, &*row_ind.begin(), col_ind.size(),
                            &*col_ind.begin(), &*locMatK.data().begin(),
                            ADD_VALUES);
      CHKERR ::MatSetValues(C, nb_rows, &*row_ind.begin(), col_ind.size(),
                            &*col_ind.begin(), &*locMatC.data().begin(),
                            ADD_VALUES);
      CHKERR ::MatSetValues(M, nb_rows, &*row_ind.begin(), col_ind.size(),
                            &*col_ind.begin(), &*locMatM.data().begin(),
                            ADD_VALUES);
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  // iterate the sides rows
  for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
    const auto sense_row = senseMap[s0];
 
    for (auto x0 = 0; x0 != indices_side_map[s0].size(); ++x0) {
 
      for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
        const auto sense_col = senseMap[s1];
 
        for (auto x1 = 0; x1 != indices_side_map[s1].size(); ++x1) {
 
          CHKERR integrate(sense_row, indices_side_map[s0][x0],
                           base_side_map[s0][x0], diff_base_side_map[s0][x0],
 
                           sense_col, indices_side_map[s1][x1],
                           base_side_map[s1][x1], diff_base_side_map[s1][x1],
 
                           *d_ptr[s0], *d_ptr[s1], *d_vis_ptr[s0],
                           *d_vis_ptr[s1]
 
          );
        }
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}