MixTransport::MixTransportElement::OpEvaluateBcOnFluxes Struct Reference

Evaluate boundary conditions on fluxes. More...

#include <users_modules/tutorials/cor-0to1/src/MixTransportElement.hpp>

Inheritance diagram for MixTransport::MixTransportElement::OpEvaluateBcOnFluxes:

Collaboration diagram for MixTransport::MixTransportElement::OpEvaluateBcOnFluxes:

Public Member Functions
	OpEvaluateBcOnFluxes (MixTransportElement &ctx, const std::string flux_name)
MoFEMErrorCode	doWork (int side, EntityType type, DataForcesAndSourcesCore::EntData &data)

Public Attributes
MixTransportElement &	cTx
MatrixDouble	NN
VectorDouble	Nf
FTensor::Index< 'i', 3 >	i

Detailed Description

Evaluate boundary conditions on fluxes.

Note that Neumann boundary conditions here are essential. So it is opposite what you find in displacement finite element method.

Here we have to solve for degrees of freedom on boundary such base functions approximate flux.

Definition at line 1126 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

OpEvaluateBcOnFluxes()

MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::OpEvaluateBcOnFluxes	(	MixTransportElement &	ctx,
		const std::string	flux_name
	)

Definition at line 1129 of file MixTransportElement.hpp.

        : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
              flux_name, UserDataOperator::OPROW),
          cTx(ctx) {}

Member Function Documentation

doWork()

MoFEMErrorCode MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::doWork	(	int	side,
		EntityType	type,
		DataForcesAndSourcesCore::EntData &	data
	)

Definition at line 1138 of file MixTransportElement.hpp.

                                                                 {
      MoFEMFunctionBegin;
      if (data.getFieldData().size() == 0)
        MoFEMFunctionReturnHot(0);
      EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
      int nb_dofs = data.getFieldData().size();
      int nb_gauss_pts = data.getN().size1();
      if (3 * nb_dofs != static_cast<int>(data.getN().size2())) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "wrong number of dofs");
      }
      NN.resize(nb_dofs, nb_dofs);
      Nf.resize(nb_dofs);
      NN.clear();
      Nf.clear();
 
      // Get normal vector. Note that when higher order geometry is set, then
      // face element could be curved, i.e. normal can be different at each
      // integration point.
      double *normal_ptr;
      if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
        // HO geometry
        normal_ptr = &getNormalsAtGaussPts(0)[0];
      } else {
        // Linear geometry, i.e. constant normal on face
        normal_ptr = &getNormal()[0];
      }
      // set tensor from pointer
      FTensor::Tensor1<const double *, 3> t_normal(normal_ptr, &normal_ptr[1],
                                                   &normal_ptr[2], 3);
      // get base functions
      auto t_n_hdiv_row = data.getFTensor1N<3>();
 
      double nrm2 = 0;
 
      // loop over integration points
      for (int gg = 0; gg < nb_gauss_pts; gg++) {
 
        // get integration point coordinates
        const double x = getCoordsAtGaussPts()(gg, 0);
        const double y = getCoordsAtGaussPts()(gg, 1);
        const double z = getCoordsAtGaussPts()(gg, 2);
 
        // get flux on fece for given element handle and coordinates
        double flux;
        CHKERR cTx.getBcOnFluxes(fe_ent, x, y, z, flux);
        ;
        // get weight for integration rule
        double w = getGaussPts()(2, gg);
        if (gg == 0) {
          nrm2 = sqrt(t_normal(i) * t_normal(i));
        }
 
        // set tensor of rank 0 to matrix NN elements
        // loop over base functions on rows and columns
        for (int ll = 0; ll != nb_dofs; ll++) {
          // get column on shape functions
          FTensor::Tensor1<const double *, 3> t_n_hdiv_col(
              &data.getVectorN<3>(gg)(0, HVEC0),
              &data.getVectorN<3>(gg)(0, HVEC1),
              &data.getVectorN<3>(gg)(0, HVEC2), 3);
          for (int kk = 0; kk <= ll; kk++) {
            NN(ll, kk) += w * t_n_hdiv_row(i) * t_n_hdiv_col(i);
            ++t_n_hdiv_col;
          }
          // right hand side
          Nf[ll] += w * t_n_hdiv_row(i) * t_normal(i) * flux / nrm2;
          ++t_n_hdiv_row;
        }
 
        // If HO geometry increment t_normal to next integration point
        if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
          ++t_normal;
          nrm2 = sqrt(t_normal(i) * t_normal(i));
        }
      }
      // get global dofs indices on element
      cTx.bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
      // factor matrix
      cholesky_decompose(NN);
      // solve local problem
      cholesky_solve(NN, Nf, ublas::lower());
 
      // set solution to vector
      CHKERR VecSetOption(cTx.D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
      CHKERR VecSetValues(cTx.D0, nb_dofs, &*data.getIndices().begin(),
                          &*Nf.begin(), INSERT_VALUES);
      for (int dd = 0; dd != nb_dofs; ++dd)
        data.getFieldDofs()[dd]->getFieldData() = Nf[dd];
 
      MoFEMFunctionReturn(0);
    }

Member Data Documentation

cTx

MixTransportElement& MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::cTx

Definition at line 1128 of file MixTransportElement.hpp.

i

FTensor::Index<'i', 3> MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::i

Definition at line 1136 of file MixTransportElement.hpp.

Nf

VectorDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::Nf

Definition at line 1135 of file MixTransportElement.hpp.

NN

MatrixDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::NN

Definition at line 1134 of file MixTransportElement.hpp.

---

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

• users_modules/tutorials/cor-0to1/src/MixTransportElement.hpp