20 double &tangent_stiffness,
21 boost::shared_ptr<Range> &ents,
25 auto cubit_meshset_ptr =
26 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(ms_id,
28 std::vector<double> block_data;
29 CHKERR cubit_meshset_ptr->getAttributes(block_data);
30 if (block_data.size() != 2) {
32 "Expected that block has two attribute");
34 normal_stiffness = block_data[0];
35 tangent_stiffness = block_data[1];
39 <<
"Normal stiffness " << normal_stiffness;
41 <<
"Tangent stiffness " << tangent_stiffness;
43 ents = boost::make_shared<Range>();
45 m_field.
get_moab().get_entities_by_handle(cubit_meshset_ptr->meshset,
58template <
int FIELD_DIM, AssemblyType A,
typename EleOp>
65 boost::shared_ptr<MatrixDouble> u_ptr,
double scale);
71 boost::shared_ptr<MatrixDouble>
uPtr;
72 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &data);
75template <
int FIELD_DIM, AssemblyType A,
typename EleOp>
82 std::string row_field_name, std::string col_field_name);
87 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &row_data,
88 EntitiesFieldData::EntData &col_data);
91template <
int FIELD_DIM, AssemblyType A,
typename EleOp>
92OpFluxRhsImpl<ElasticExample::SpringBcType<BLOCKSET>, 1,
FIELD_DIM,
A, GAUSS,
95 boost::shared_ptr<MatrixDouble> u_ptr,
99 this->normalStiffness, this->tangentStiffness,
100 this->entsPtr, m_field, ms_id),
101 "Can not read spring stiffness data from blockset");
104template <
int FIELD_DIM, AssemblyType A,
typename EleOp>
107 GAUSS, EleOp>::iNtegrate(EntitiesFieldData::EntData
113 const double vol = OpBase::getMeasure();
115 auto t_w = OpBase::getFTensor0IntegrationWeight();
117 auto t_row_base = row_data.getFTensor0N();
120 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
123 auto t_u = getFTensor1FromMat<FIELD_DIM>(*uPtr);
127 const double alpha = t_w * vol * rhsScale;
129 auto l2_norm = t_normal(
i) * t_normal(
i);
133 t_P(
i,
j) = t_normal(
i) * t_normal(
j) / l2_norm;
140 t_D(
i,
j) = normalStiffness * t_P(
i,
j) + tangentStiffness * t_Q(
i,
j);
144 t_reaction(
i) = t_D(
i,
j) * t_u(
j);
147 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(
OpBase::locF);
150 t_nf(
i) += (alpha * t_row_base) * t_reaction(
i);
165template <
int FIELD_DIM, AssemblyType A,
typename EleOp>
166OpFluxLhsImpl<ElasticExample::SpringBcType<BLOCKSET>, 1,
FIELD_DIM,
A,
GAUSS,
168 std::string row_field_name,
169 std::string col_field_name)
170 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL) {
172 this->normalStiffness, this->tangentStiffness,
173 this->entsPtr, m_field, ms_id),
174 "Can not read spring stiffness data from blockset");
177template <
int FIELD_DIM, AssemblyType A,
typename EleOp>
181 EleOp>::iNtegrate(EntitiesFieldData::EntData &row_data,
182 EntitiesFieldData::EntData &col_data) {
189 const double vol = OpBase::getMeasure();
191 auto t_w = OpBase::getFTensor0IntegrationWeight();
193 auto t_row_base = row_data.getFTensor0N();
196 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
201 const double alpha = t_w * vol;
203 auto l2_norm = t_normal(
i) * t_normal(
i);
207 t_P(
i,
j) = t_normal(
i) * t_normal(
j) / l2_norm;
214 t_D(
i,
j) = normalStiffness * t_P(
i,
j) + tangentStiffness * t_Q(
i,
j);
220 auto t_col_base = col_data.getFTensor0N(gg, 0);
222 auto t_m = OpBase::template getLocMat<FIELD_DIM>(
FIELD_DIM * rr);
226 t_m(
i,
j) += t_D(
i,
j) * (alpha * t_row_base * t_col_base);
243struct AddFluxToRhsPipelineImpl<
253 static MoFEMErrorCode
add(
255 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
257 boost::shared_ptr<MatrixDouble> u_ptr,
double scale,
258 std::string block_name, Sev sev
263 using OP = OpFluxRhsImpl<ElasticExample::SpringBcType<BLOCKSET>,
BASE_DIM,
266 auto add_op = [&](
auto &&meshset_vec_ptr) {
267 for (
auto m : meshset_vec_ptr) {
279 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
282 (boost::format(
"%s(.*)") % block_name).str()
291 "Handling of bc type not implemented");
300struct AddFluxToLhsPipelineImpl<
310 static MoFEMErrorCode
add(
312 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
314 std::string col_field_name, std::string block_name, Sev sev
319 using OP = OpFluxLhsImpl<ElasticExample::SpringBcType<BLOCKSET>,
BASE_DIM,
322 auto add_op = [&](
auto &&meshset_vec_ptr) {
323 for (
auto m : meshset_vec_ptr) {
326 new OP(m_field,
m->getMeshsetId(), row_field_name, col_field_name));
335 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
338 (boost::format(
"%s(.*)") % block_name).str()
347 "Handling of bc type not implemented");
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
CubitBC
Types of sets and boundary conditions.
@ 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< 'm', SPACE_DIM > m
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
constexpr IntegrationType I
constexpr auto field_name
static MoFEMErrorCode add(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, MoFEM::Interface &m_field, std::string row_field_name, std::string col_field_name, std::string block_name, Sev sev)
AddFluxToLhsPipelineImpl()=delete
AddFluxToRhsPipelineImpl()=delete
static MoFEMErrorCode add(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, MoFEM::Interface &m_field, std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, double scale, std::string block_name, Sev sev)
static MoFEMErrorCode getStiffness(double &normal_stiffness, double &tangent_stiffness, boost::shared_ptr< Range > &ents, MoFEM::Interface &m_field, int ms_id)
GetSpringStiffness()=delete
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< MatrixDouble > uPtr