10#ifndef __FORMSBROKENSPACECONSTRAINTIMPL_HPP__
11#define __FORMSBROKENSPACECONSTRAINTIMPL_HPP__
45 prev_side_fe_ptr->numeredEntFiniteElementPtr->getFEUId();
58 prev_side_fe_ptr->numeredEntFiniteElementPtr;
63 "sideFEName is different");
92 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data);
101template <
typename OpBase>
104 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data)
106 brokenBaseSideData(broken_base_side_data) {}
108template <
typename OpBase>
117 brokenBaseSideData->resize(OP::getLoopSize());
119 const auto n_in_the_loop = OP::getNinTheLoop();
120 const auto face_sense = OP::getSkeletonSense();
123 if (face_sense != -1 && face_sense != 1)
127 auto set_data = [&](
auto &side_data) {
128 side_data.getSide() = row_side;
129 side_data.getType() = row_type;
130 side_data.getSense() = face_sense;
131 side_data.getData().sEnse = row_data.
sEnse;
132 side_data.getData().sPace = row_data.
sPace;
133 side_data.getData().bAse = row_data.
bAse;
134 side_data.getData().iNdices = row_data.
iNdices;
135 side_data.getData().localIndices = row_data.
localIndices;
136 side_data.getData().dOfs = row_data.
dOfs;
138 side_data.getData().fieldData = row_data.
fieldData;
141 auto set_base = [&](
auto &side_data) {
142 auto base = side_data.getData().getBase();
143 for (
auto dd = 0; dd != BaseDerivatives::LastDerivative; ++dd) {
144 side_data.getData().baseFunctionsAndBaseDerivatives[dd][base].reset();
146 side_data.getData().baseFunctionsAndBaseDerivatives[dd][base] =
147 boost::make_shared<MatrixDouble>(*base_ptr);
152 set_data((*brokenBaseSideData)[n_in_the_loop]);
153 set_base((*brokenBaseSideData)[n_in_the_loop]);
163 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
164 boost::shared_ptr<MatrixDouble> flux_ptr);
174template <
typename OpBase>
176 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
177 boost::shared_ptr<MatrixDouble> flux_ptr)
178 :
OP(
NOSPACE,
OP::OPSPACE), brokenBaseSideData(broken_base_side_data),
181template <
typename OpBase>
185 auto swap_flux = [&](
auto &side_data) { side_data.getFlux().swap(*fluxPtr); };
186 swap_flux((*brokenBaseSideData)[OP::getNinTheLoop()]);
195 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
196 boost::shared_ptr<Range> ents_ptr =
nullptr)
202 const std::string row_field,
203 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
204 const bool assmb_transpose,
const bool only_transpose,
205 boost::shared_ptr<Range> ents_ptr =
nullptr)
206 :
OP(row_field, row_field,
OP::OPROW, ents_ptr),
221template <
typename OpBase>
228 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
233 if (!brokenBaseSideData) {
238 auto do_work_rhs = [
this](
int row_side, EntityType row_type,
246 OP::nbIntegrationPts = OP::getGaussPts().size2();
248 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
250 OP::locF.resize(OP::nbRows,
false);
253 CHKERR this->iNtegrate(row_data);
256 CHKERR this->aSsemble(row_data);
260 auto do_work_lhs = [
this](
int row_side,
int col_side, EntityType row_type,
266 auto check_if_assemble_transpose = [&] {
268 if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
272 }
else if (OP::assembleTranspose) {
278 OP::rowSide = row_side;
279 OP::rowType = row_type;
280 OP::colSide = col_side;
281 OP::colType = col_type;
282 OP::nbCols = col_data.getIndices().size();
283 OP::locMat.resize(OP::nbRows, OP::nbCols,
false);
285 CHKERR this->iNtegrate(row_data, col_data);
287 CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
291 switch (OP::opType) {
297 OP::nbIntegrationPts = OP::getGaussPts().size2();
298 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
303 for (
auto &bd : *brokenBaseSideData) {
306 if (bd.getData().getSpace() !=
HDIV && bd.getData().getSpace() !=
HCURL) {
308 (std::string(
"Expect Hdiv or Hcurl space but received ") +
312 if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
314 "base functions not set");
321 row_side, bd.getSide(),
324 row_type, bd.getType(),
327 row_data, bd.getData(),
337 for (
auto &bd : *brokenBaseSideData) {
338 CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData(),
344 (std::string(
"wrong op type ") +
352template <
int FIELD_DIM, IntegrationType I,
typename OpBrokenBase>
355template <
int FIELD_DIM,
typename OpBrokenBase>
362 const std::string row_field,
363 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
364 boost::shared_ptr<double> beta_ptr,
const bool assmb_transpose,
365 const bool only_transpose, boost::shared_ptr<Range> ents_ptr =
nullptr)
366 :
OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
368 scalarBetaPtr(beta_ptr) {}
371 const std::string row_field,
372 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
373 double beta,
const bool assmb_transpose,
const bool only_transpose,
374 boost::shared_ptr<Range> ents_ptr =
nullptr)
376 boost::make_shared<
double>(beta),
377 assmb_transpose, only_transpose, ents_ptr) {}
385template <
int FIELD_DIM,
typename OpBase>
391 auto nb_row_dofs = row_data.
getIndices().size();
392 auto nb_col_dofs = col_data.
getIndices().size();
393 if (!nb_row_dofs || !nb_col_dofs)
399 auto t_w = this->getFTensor0IntegrationWeight();
400 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
401 size_t nb_base_functions = col_data.
getN().size2() / 3;
406 "number of dofs not divisible by field dimension");
411 "number of dofs exceeds number of base functions");
420 for (; cc != nb_col_dofs /
FIELD_DIM; cc++) {
422 for (
auto rr = 0; rr != nb_row_dofs /
FIELD_DIM; ++rr) {
424 (0.5 * t_w * t_row_base) * (t_normal(
J) * t_col_base(
J));
431 for (; cc < nb_base_functions; ++cc)
438 for (
auto rr = 0; rr != nb_row_dofs /
FIELD_DIM; ++rr) {
439 for (
auto cc = 0; cc != nb_col_dofs /
FIELD_DIM; ++cc) {
448 OP::locMat *= *scalarBetaPtr;
453template <
int FIELD_DIM, IntegrationType I,
typename OpBase>
456template <
int FIELD_DIM, IntegrationType I,
typename OpBase>
459template <
int FIELD_DIM,
typename OpBrokenBase>
466 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
467 boost::shared_ptr<MatrixDouble> lagrange_ptr,
468 boost::shared_ptr<double> beta_ptr,
469 boost::shared_ptr<Range> ents_ptr =
nullptr)
470 :
OP(broken_base_side_data, ents_ptr), scalarBetaPtr(beta_ptr),
471 lagrangePtr(lagrange_ptr) {}
474 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
475 boost::shared_ptr<MatrixDouble> lagrange_ptr,
double beta,
476 boost::shared_ptr<Range> ents_ptr =
nullptr)
478 boost::make_shared<
double>(beta),
488template <
int FIELD_DIM,
typename OpBase>
497 auto t_w = this->getFTensor0IntegrationWeight();
498 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
502 auto nb_base_functions = row_data.
getN().size2() / 3;
508 t_vec(
i) += (0.5 * t_w * (t_row_base(
J) * t_normal(
J))) * t_lagrange(
i);
512 for (; rr < nb_base_functions; ++rr)
520 OP::locF *= *scalarBetaPtr;
525template <
int FIELD_DIM,
typename OpBase>
532 const std::string row_field,
533 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_side_data_ptr,
534 boost::shared_ptr<double> beta_ptr,
535 boost::shared_ptr<Range> ents_ptr =
nullptr)
536 :
OpBase(row_field, row_field,
OpBase::OPROW, ents_ptr),
537 brokenSideDataPtr(broken_side_data_ptr), scalarBetaPtr(beta_ptr) {}
540 const std::string row_field,
541 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_side_data_ptr,
542 double beta, boost::shared_ptr<Range> ents_ptr =
nullptr)
544 boost::make_shared<
double>(beta),
553template <
int FIELD_DIM,
typename OpBase>
562 OP::locF.resize(row_data.
getIndices().size(),
false);
565 for (
auto &bd : *brokenSideDataPtr) {
566 auto t_w = this->getFTensor0IntegrationWeight();
567 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
570 auto nb_base_functions = row_data.
getN().size2() / 3;
571 auto sense = bd.getSense();
577 (sense * 0.5 * t_w) * t_row_base * t_normal(
J) * t_flux(
i,
J);
581 for (; rr < nb_base_functions; ++rr)
590 OP::locF *= *scalarBetaPtr;
#define FTENSOR_INDEX(DIM, I)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
static const char *const FieldSpaceNames[]
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
EntitiesFieldData::EntData entData
Data on single entity (This is passed as argument to DataOperator::doWork)
FieldSpace sPace
Entity space.
int sEnse
Entity sense (orientation)
VectorDouble fieldData
Field data on entity.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
VectorInt localIndices
Local indices on entity.
VectorInt iNdices
Global indices on entity.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
virtual int getSense() const
get entity sense, need to calculate base functions with conforming approximation fields
VectorFieldEntities fieldEntities
Field entities.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
std::array< std::array< boost::shared_ptr< MatrixDouble >, LASTBASE >, LastDerivative > baseFunctionsAndBaseDerivatives
VectorDofs dOfs
DoFs on entity.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FieldApproximationBase bAse
Field approximation base.
static const char *const OpTypeNames[]
ForcesAndSourcesCore * getPtrFE() const
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
ForcesAndSourcesCore * getSidePtrFE() const
boost::shared_ptr< Range > entsPtr
Entities on which element is run.
int nbIntegrationPts
number of integration points
boost::shared_ptr< std::vector< BrokenBaseSideData > > brokenBaseSideData
OpBrokenBaseImpl(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, const bool assmb_transpose, const bool only_transpose, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpBrokenBaseImpl(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< Range > ents_ptr=nullptr)
Operator for broken loop side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< double > scalarBetaPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
OpBrokenSpaceConstrainDFluxImpl(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< MatrixDouble > lagrange_ptr, boost::shared_ptr< double > beta_ptr, boost::shared_ptr< Range > ents_ptr=nullptr)
OpBrokenSpaceConstrainDFluxImpl(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< MatrixDouble > lagrange_ptr, double beta, boost::shared_ptr< Range > ents_ptr=nullptr)
boost::shared_ptr< MatrixDouble > lagrangePtr
OpBrokenSpaceConstrainDHybridImpl(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_side_data_ptr, boost::shared_ptr< double > beta_ptr, boost::shared_ptr< Range > ents_ptr=nullptr)
OpBrokenSpaceConstrainDHybridImpl(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_side_data_ptr, double beta, boost::shared_ptr< Range > ents_ptr=nullptr)
boost::shared_ptr< double > scalarBetaPtr
boost::shared_ptr< std::vector< BrokenBaseSideData > > brokenSideDataPtr
boost::shared_ptr< double > scalarBetaPtr
OpBrokenSpaceConstrainImpl(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< double > beta_ptr, const bool assmb_transpose, const bool only_transpose, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpBrokenSpaceConstrainImpl(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, double beta, const bool assmb_transpose, const bool only_transpose, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
boost::shared_ptr< std::vector< BrokenBaseSideData > > brokenBaseSideData
OpGetBrokenBaseSideData(const std::string field_name, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data)
Element used to execute operators on side of the element.
const std::string sideFEName
boost::shared_ptr< E > sideFEPtr
OpLoopSide(MoFEM::Interface &m_field, const std::string fe_name, const int side_dim, boost::shared_ptr< Range > fe_range, const LogManager::SeverityLevel sev=Sev::noisy, boost::shared_ptr< AdjCache > adj_cache=nullptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpSetFlux(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< MatrixDouble > flux_ptr)
boost::shared_ptr< std::vector< BrokenBaseSideData > > brokenBaseSideData
boost::shared_ptr< MatrixDouble > fluxPtr