#include <users_modules/basic_finite_elements/src/ConvectiveMassElement.hpp>

Inheritance diagram for ConvectiveMassElement::OpMassJacobian:

Collaboration diagram for ConvectiveMassElement::OpMassJacobian:

Public Member Functions
	OpMassJacobian (const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)

MoFEMErrorCode	doWork (int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)

Public Attributes
BlockData &	dAta

CommonData &	commonData

int	tAg

bool	jAcobian

bool &	lInear

bool	fieldDisp

boost::ptr_vector< MethodForForceScaling > &	methodsOp

FTensor::Index< 'i', 3 >	i

FTensor::Index< 'j', 3 >	j

FTensor::Index< 'k', 3 >	k

VectorBoundedArray< adouble, 3 >	a

VectorBoundedArray< adouble, 3 >	dot_W

VectorBoundedArray< adouble, 3 >	dp_dt

VectorBoundedArray< adouble, 3 >	a_res

MatrixBoundedArray< adouble, 9 >	h

MatrixBoundedArray< adouble, 9 >	H

MatrixBoundedArray< adouble, 9 >	invH

MatrixBoundedArray< adouble, 9 >	F

MatrixBoundedArray< adouble, 9 >	g

MatrixBoundedArray< adouble, 9 >	G

std::vector< double >	active

Detailed Description

Definition at line 177 of file ConvectiveMassElement.hpp.

Constructor & Destructor Documentation

◆ OpMassJacobian()

ConvectiveMassElement::OpMassJacobian::OpMassJacobian	(	const std::string	field_name,
		BlockData &	data,
		CommonData &	common_data,
		boost::ptr_vector< MethodForForceScaling > &	methods_op,
		int	tag,
		bool	linear = `false`
	)

Definition at line 154 of file ConvectiveMassElement.cpp.

    : VolumeElementForcesAndSourcesCore::UserDataOperator(
          field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
      dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
      lInear(commonData.lInear), fieldDisp(false), methodsOp(methods_op) {}

Member Function Documentation

◆ doWork()

MoFEMErrorCode ConvectiveMassElement::OpMassJacobian::doWork	(	int	row_side,
		EntityType	row_type,
		EntitiesFieldData::EntData &	row_data
	)

Definition at line 163 of file ConvectiveMassElement.cpp.

                                                                           {
  MoFEMFunctionBeginHot;
 
  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
      dAta.tEts.end()) {
    MoFEMFunctionReturnHot(0);
  }
 
  // do it only once, no need to repeat this for edges,faces or tets
  if (row_type != MBVERTEX)
    MoFEMFunctionReturnHot(0);
 
  int nb_dofs = row_data.getIndices().size();
  if (nb_dofs == 0)
    MoFEMFunctionReturnHot(0);
 
  {
 
    if (a.size() != 3) {
      a.resize(3, false);
      dot_W.resize(3, false);
      a_res.resize(3, false);
      g.resize(3, 3, false);
      G.resize(3, 3, false);
      h.resize(3, 3, false);
      H.resize(3, 3, false);
      invH.resize(3, 3, false);
      F.resize(3, 3, false);
    }
 
    std::fill(dot_W.begin(), dot_W.end(), 0);
    std::fill(H.data().begin(), H.data().end(), 0);
    std::fill(invH.data().begin(), invH.data().end(), 0);
    for (int ii = 0; ii != 3; ii++) {
      H(ii, ii) = 1;
      invH(ii, ii) = 1;
    }
 
    int nb_gauss_pts = row_data.getN().size1();
    commonData.valMass.resize(nb_gauss_pts);
    commonData.jacMassRowPtr.resize(nb_gauss_pts);
    commonData.jacMass.resize(nb_gauss_pts);
 
    const std::vector<VectorDouble> &dot_spacial_vel =
        commonData.dataAtGaussPts["DOT_" + commonData.spatialVelocities];
 
    const std::vector<MatrixDouble> &spatial_positions_grad =
        commonData.gradAtGaussPts[commonData.spatialPositions];
 
    const std::vector<MatrixDouble> &spatial_velocities_grad =
        commonData.gradAtGaussPts[commonData.spatialVelocities];
 
    const std::vector<VectorDouble> &meshpos_vel =
        commonData.dataAtGaussPts["DOT_" + commonData.meshPositions];
 
    const std::vector<MatrixDouble> &mesh_positions_gradient =
        commonData.gradAtGaussPts[commonData.meshPositions];
 
    int nb_active_vars = 0;
    for (int gg = 0; gg < nb_gauss_pts; gg++) {
 
      if (gg == 0) {
 
        trace_on(tAg);
 
        for (int nn1 = 0; nn1 < 3; nn1++) { // 0
          // commonData.dataAtGaussPts["DOT_"+commonData.spatialVelocities]
          a[nn1] <<= dot_spacial_vel[gg][nn1];
          nb_active_vars++;
        }
        for (int nn1 = 0; nn1 < 3; nn1++) { // 3
          for (int nn2 = 0; nn2 < 3; nn2++) {
            // commonData.gradAtGaussPts[commonData.spatialPositions][gg]
            h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
            if (fieldDisp) {
              if (nn1 == nn2) {
                h(nn1, nn2) += 1;
              }
            }
            nb_active_vars++;
          }
        }
        if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions]
                .size() > 0) {
          for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
            for (int nn2 = 0; nn2 < 3; nn2++) {
              // commonData.gradAtGaussPts[commonData.spatialVelocities]
              g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
              nb_active_vars++;
            }
          }
          for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
            // commonData.dataAtGaussPts["DOT_"+commonData.meshPositions]
            dot_W(nn1) <<= meshpos_vel[gg][nn1];
            nb_active_vars++;
          }
          for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
            for (int nn2 = 0; nn2 < 3; nn2++) {
              // commonData.gradAtGaussPts[commonData.meshPositions][gg]
              H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
              nb_active_vars++;
            }
          }
        }
 
        auto a0 = dAta.a0;
        CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, a0);
 
        auto t_a_res =
            FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
        auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
        auto t_a0 = FTensor::Tensor1<double *, 3>{&a0[0], &a0[1], &a0[2]};
        auto t_dotW =
            FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
        auto t_g = getFTensor2FromArray3by3(g, FTensor::Number<0>(), 0);
        auto t_G = getFTensor2FromArray3by3(G, FTensor::Number<0>(), 0);
        auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
        auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
        auto t_h = getFTensor2FromArray3by3(h, FTensor::Number<0>(), 0);
 
        const double rho0 = dAta.rho0;
 
        adouble detH = determinantTensor3by3(H);
        CHKERR invertTensor3by3(H, detH, invH);
 
        t_G(i, j) = t_g(i, k) * t_invH(k, j);
        t_a_res(i) = t_a(i) - t_a0(i) + t_G(i, j) * t_dotW(j);
 
        // FIXME: there is error somewhere for nonlinear case
        //  test dam example with -is_linear 0
        if (!lInear) {
 
          t_F(i, j) = t_h(i, k) * t_invH(k, j);
          t_a_res(i) *= rho0 * detH;
          t_a_res(i) *= determinantTensor3by3(t_F);
 
        } else {
 
          t_a_res(i) *= rho0 * detH;
        }
 
        // dependant
        VectorDouble &res = commonData.valMass[gg];
        res.resize(3);
        for (int rr = 0; rr < 3; rr++) {
          a_res[rr] >>= res[rr];
        }
 
        trace_off();
      }
 
      active.resize(nb_active_vars);
      int aa = 0;
      for (int nn1 = 0; nn1 < 3; nn1++) { // 0
        active[aa++] = dot_spacial_vel[gg][nn1];
      }
      for (int nn1 = 0; nn1 < 3; nn1++) { // 3
        for (int nn2 = 0; nn2 < 3; nn2++) {
          if (fieldDisp && nn1 == nn2) {
            active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
          } else {
            active[aa++] = spatial_positions_grad[gg](nn1, nn2);
          }
        }
      }
      if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
          0) {
        for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
          for (int nn2 = 0; nn2 < 3; nn2++) {
            active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
          }
        }
        for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
          active[aa++] = meshpos_vel[gg][nn1];
        }
        for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
          for (int nn2 = 0; nn2 < 3; nn2++) {
            active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
          }
        }
      }
 
      if (!jAcobian) {
        VectorDouble &res = commonData.valMass[gg];
        if (gg > 0) {
          res.resize(3);
          int r;
          r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
          if (r != 3) { // function is locally analytic
            SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                     "ADOL-C function evaluation with error r = %d", r);
          }
        }
        double val = getVolume() * getGaussPts()(3, gg);
        res *= val;
        // cout << "my res " << res << endl;
      } else {
        commonData.jacMassRowPtr[gg].resize(3);
        commonData.jacMass[gg].resize(3, nb_active_vars);
        for (int nn1 = 0; nn1 < 3; nn1++) {
          (commonData.jacMassRowPtr[gg])[nn1] =
              &(commonData.jacMass[gg](nn1, 0));
        }
        int r;
        r = jacobian(tAg, 3, nb_active_vars, &active[0],
                     &(commonData.jacMassRowPtr[gg])[0]);
        if (r != 3) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                  "ADOL-C function evaluation with error");
        }
        double val = getVolume() * getGaussPts()(3, gg);
        commonData.jacMass[gg] *= val;
      }
    }
  }
 
  MoFEMFunctionReturnHot(0);
}

Member Data Documentation

◆ a

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::a

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ a_res

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::a_res

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ active

std::vector<double> ConvectiveMassElement::OpMassJacobian::active

Definition at line 201 of file ConvectiveMassElement.hpp.

◆ commonData

CommonData& ConvectiveMassElement::OpMassJacobian::commonData

Definition at line 182 of file ConvectiveMassElement.hpp.

◆ dAta

BlockData& ConvectiveMassElement::OpMassJacobian::dAta

Definition at line 181 of file ConvectiveMassElement.hpp.

◆ dot_W

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::dot_W

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ dp_dt

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::dp_dt

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ F

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::F

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ fieldDisp

bool ConvectiveMassElement::OpMassJacobian::fieldDisp

Definition at line 186 of file ConvectiveMassElement.hpp.

◆ g

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::g

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ G

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::G

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ h

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::h

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ H

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::H

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ i

FTensor::Index<'i', 3> ConvectiveMassElement::OpMassJacobian::i

Definition at line 195 of file ConvectiveMassElement.hpp.

◆ invH

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::invH

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ j

FTensor::Index<'j', 3> ConvectiveMassElement::OpMassJacobian::j

Definition at line 196 of file ConvectiveMassElement.hpp.

◆ jAcobian

bool ConvectiveMassElement::OpMassJacobian::jAcobian

Definition at line 184 of file ConvectiveMassElement.hpp.

◆ k

FTensor::Index<'k', 3> ConvectiveMassElement::OpMassJacobian::k

Definition at line 197 of file ConvectiveMassElement.hpp.

◆ lInear

bool& ConvectiveMassElement::OpMassJacobian::lInear

Definition at line 185 of file ConvectiveMassElement.hpp.

◆ methodsOp

boost::ptr_vector<MethodForForceScaling>& ConvectiveMassElement::OpMassJacobian::methodsOp

Definition at line 188 of file ConvectiveMassElement.hpp.

◆ tAg

int ConvectiveMassElement::OpMassJacobian::tAg

Definition at line 183 of file ConvectiveMassElement.hpp.

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

users_modules/basic_finite_elements/src/ConvectiveMassElement.hpp
users_modules/basic_finite_elements/src/impl/ConvectiveMassElement.cpp

Public Member Functions

Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ OpMassJacobian()

Member Function Documentation

◆ doWork()

Member Data Documentation

◆ a

◆ a_res

◆ active

◆ commonData

◆ dAta

◆ dot_W

◆ dp_dt

◆ F

◆ fieldDisp

◆ g

◆ G

◆ h

◆ H

◆ i

◆ invH

◆ j

◆ jAcobian

◆ k

◆ lInear

◆ methodsOp

◆ tAg