11#include <MethodForForceScaling.hpp>
13#include <boost/shared_ptr.hpp>
14#include <boost/numeric/ublas/vector_proxy.hpp>
22 int row_side,
int col_side, EntityType row_type, EntityType col_type,
33 const int rank = dof_ptr->getNbOfCoeffs();
35 int nb_row_dofs = row_data.
getIndices().size() / rank;
36 int nb_col_dofs = col_data.
getIndices().size() / rank;
38 NN.resize(nb_row_dofs, nb_col_dofs);
41 unsigned int nb_gauss_pts = row_data.
getN().size1();
42 for (
unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
44 double w = getGaussPts()(2, gg);
46 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
49 cblas_dger(CblasRowMajor, nb_row_dofs, nb_col_dofs, w,
50 &row_data.
getN()(gg, 0), 1, &col_data.
getN()(gg, 0), 1,
51 &*NN.data().begin(), nb_col_dofs);
54 if ((row_type != col_type) || (row_side != col_side)) {
55 transNN.resize(nb_col_dofs, nb_row_dofs);
56 ublas::noalias(transNN) = trans(NN);
59 double *data = &*NN.data().begin();
60 double *trans_data = &*transNN.data().begin();
62 row_indices.resize(nb_row_dofs);
63 col_indices.resize(nb_col_dofs);
65 for (
int rr = 0; rr < rank; rr++) {
67 if ((row_data.
getIndices().size() % rank) != 0) {
71 if ((col_data.
getIndices().size() % rank) != 0) {
82 ublas::noalias(row_indices) = ublas::vector_slice<VectorInt>(
84 ublas::slice(rr, rank, row_data.
getIndices().size() / rank));
85 ublas::noalias(col_indices) = ublas::vector_slice<VectorInt>(
87 ublas::slice(rr, rank, col_data.
getIndices().size() / rank));
89 nb_rows = row_indices.size();
90 nb_cols = col_indices.size();
91 rows = &*row_indices.data().begin();
92 cols = &*col_indices.data().begin();
102 if (nb_rows != NN.size1()) {
105 if (nb_cols != NN.size2()) {
111 if ((row_type != col_type) || (row_side != col_side)) {
112 if (nb_rows != transNN.size2()) {
114 "data inconsistency");
116 if (nb_cols != transNN.size1()) {
118 "data inconsistency");
121 trans_data, ADD_VALUES);
129 const std::string &field, Mat A,
134 const std::string &field)
139 if (mapZeroRows.empty()) {
142 "Need to initialized from AnalyticalDirichletBC::solveProblem");
144 CHKERR iNitialize(*trisPtr);
151 ParallelComm *pcomm =
154 CHKERR mField.get_moab().get_connectivity(tris, ents,
true);
155 CHKERR mField.get_moab().get_adjacencies(tris, 1,
false, ents,
156 moab::Interface::UNION);
159 const auto bit_number = mField.get_field_bit_number(fieldName);
161 for (
auto eit = ents.pair_begin(); eit != ents.pair_end(); eit++) {
162 const auto f = eit->first;
163 const auto s = eit->second;
165 auto &dofs = *problemPtr->getNumeredRowDofsPtr();
170 for (; dit != hi_dit; ++dit) {
171 if ((*dit)->getPart() == mField.get_comm_rank()) {
172 mapZeroRows[(*dit)->getPetscGlobalDofIdx()] = (*dit)->getFieldData();
177 dofsIndices.resize(mapZeroRows.size());
178 dofsValues.resize(mapZeroRows.size());
180 std::map<DofIdx, FieldData>::iterator mit = mapZeroRows.begin();
181 for (; mit != mapZeroRows.end(); mit++, ii++) {
182 dofsIndices[ii] = mit->first;
183 dofsValues[ii] = mit->second;
193 string field,
Range &tris,
194 string nodals_positions) {
214 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(problem, &A);
223 string problem,
string fe,
235 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
236 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
239 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
240 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
243 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
244 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
247 problem,
ROW,
D, INSERT_VALUES, SCATTER_REVERSE);
258 string problem,
string fe,
263 CHKERR m_field.
get_moab().get_entities_by_type(fe_meshset, MBTRI, bc_tris);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< int > VectorInt
implementation of Data Operators for Forces and Sources
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
constexpr auto field_name
OpLhs(const std::string field_name)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MyTriFE & getLoopFeApprox()
Structure used to enforce analytical boundary conditions.
MoFEMErrorCode iNitialize()
DirichletBC(MoFEM::Interface &m_field, const std::string &field, Mat A, Vec X, Vec F)
boost::shared_ptr< Range > trisPtr
AnalyticalDirichletBC(MoFEM::Interface &m_field)
MoFEMErrorCode destroyProblem()
Destroy problem.
MoFEMErrorCode setUpProblem(MoFEM::Interface &m_field, string problem)
set problem solver and create matrices and vectors
MoFEMErrorCode setFiniteElement(MoFEM::Interface &m_field, string fe, string field, Range &tris, string nodals_positions="MESH_NODE_POSITIONS")
set finite element
MoFEMErrorCode solveProblem(MoFEM::Interface &m_field, string problem, string fe, DirichletBC &bc, Range &tris)
solve boundary problem
Set Dirichlet boundary conditions on displacements.
std::map< DofIdx, FieldData > mapZeroRows
std::vector< int > dofsIndices
std::vector< double > dofsValues
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure to get information form mofem into EntitiesFieldData
Matrix manager is used to build and partition problems.
Mat & snes_B
preconditioner of jacobian matrix
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.