v0.14.0
Loading...
Searching...
No Matches
AnalyticalDirichlet.cpp
Go to the documentation of this file.
1/** \file AnalyticalDirichlet.cpp
2
3 Enforce Dirichlet boundary condition for given analytical function,
4
5*/
6
7
8
9#include <MoFEM.hpp>
10using namespace MoFEM;
12#include <DirichletBC.hpp>
13#include <boost/shared_ptr.hpp>
14#include <boost/numeric/ublas/vector_proxy.hpp>
16
20
22 int row_side, int col_side, EntityType row_type, EntityType col_type,
26
27 if (row_data.getIndices().size() == 0)
29 if (col_data.getIndices().size() == 0)
31
32 const auto &dof_ptr = row_data.getFieldDofs()[0];
33 const int rank = dof_ptr->getNbOfCoeffs();
34
35 int nb_row_dofs = row_data.getIndices().size() / rank;
36 int nb_col_dofs = col_data.getIndices().size() / rank;
37
38 NN.resize(nb_row_dofs, nb_col_dofs);
39 NN.clear();
40
41 unsigned int nb_gauss_pts = row_data.getN().size1();
42 for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
43
44 double w = getGaussPts()(2, gg);
45 // higher order element
46 double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
47 w *= area;
48
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);
52 }
53
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);
57 }
58
59 double *data = &*NN.data().begin();
60 double *trans_data = &*transNN.data().begin();
61 VectorInt row_indices, col_indices;
62 row_indices.resize(nb_row_dofs);
63 col_indices.resize(nb_col_dofs);
64
65 for (int rr = 0; rr < rank; rr++) {
66
67 if ((row_data.getIndices().size() % rank) != 0) {
68 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
69 }
70
71 if ((col_data.getIndices().size() % rank) != 0) {
72 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
73 }
74
75 unsigned int nb_rows;
76 unsigned int nb_cols;
77 int *rows;
78 int *cols;
79
80 if (rank > 1) {
81
82 ublas::noalias(row_indices) = ublas::vector_slice<VectorInt>(
83 row_data.getIndices(),
84 ublas::slice(rr, rank, row_data.getIndices().size() / rank));
85 ublas::noalias(col_indices) = ublas::vector_slice<VectorInt>(
86 col_data.getIndices(),
87 ublas::slice(rr, rank, col_data.getIndices().size() / rank));
88
89 nb_rows = row_indices.size();
90 nb_cols = col_indices.size();
91 rows = &*row_indices.data().begin();
92 cols = &*col_indices.data().begin();
93
94 } else {
95
96 nb_rows = row_data.getIndices().size();
97 nb_cols = col_data.getIndices().size();
98 rows = &*row_data.getIndices().data().begin();
99 cols = &*col_data.getIndices().data().begin();
100 }
101
102 if (nb_rows != NN.size1()) {
103 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
104 }
105 if (nb_cols != NN.size2()) {
106 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
107 }
108
109 CHKERR MatSetValues(getFEMethod()->snes_B, nb_rows, rows, nb_cols, cols,
110 data, ADD_VALUES);
111 if ((row_type != col_type) || (row_side != col_side)) {
112 if (nb_rows != transNN.size2()) {
113 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
114 "data inconsistency");
115 }
116 if (nb_cols != transNN.size1()) {
117 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
118 "data inconsistency");
119 }
120 CHKERR MatSetValues(getFEMethod()->snes_B, nb_cols, cols, nb_rows, rows,
121 trans_data, ADD_VALUES);
122 }
123 }
124
126}
127
129 const std::string &field, Mat A,
130 Vec X, Vec F)
131 : DirichletDisplacementBc(m_field, field, A, X, F) {}
132
134 const std::string &field)
135 : DirichletDisplacementBc(m_field, field) {}
136
139 if (mapZeroRows.empty()) {
140 if (!trisPtr) {
141 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
142 "Need to initialized from AnalyticalDirichletBC::solveProblem");
143 }
144 CHKERR iNitialize(*trisPtr);
145 }
147}
148
151 ParallelComm *pcomm =
152 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
153 Range ents;
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);
157 ents.merge(tris);
158
159 const auto bit_number = mField.get_field_bit_number(fieldName);
160
161 for (auto eit = ents.pair_begin(); eit != ents.pair_end(); eit++) {
162 const auto f = eit->first;
163 const auto s = eit->second;
164
165 auto &dofs = *problemPtr->getNumeredRowDofsPtr();
166 auto dit = dofs.get<Unique_mi_tag>().lower_bound(
167 DofEntity::getLoFieldEntityUId(bit_number, f));
168 auto hi_dit = dofs.get<Unique_mi_tag>().upper_bound(
169 DofEntity::getHiFieldEntityUId(bit_number, s));
170 for (; dit != hi_dit; ++dit) {
171 if ((*dit)->getPart() == mField.get_comm_rank()) {
172 mapZeroRows[(*dit)->getPetscGlobalDofIdx()] = (*dit)->getFieldData();
173 }
174 }
175 }
176
177 dofsIndices.resize(mapZeroRows.size());
178 dofsValues.resize(mapZeroRows.size());
179 int ii = 0;
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;
184 }
186}
187
189 : approxField(m_field){};
190
193 string field, Range &tris,
194 string nodals_positions) {
196
197 CHKERR m_field.add_finite_element(fe, MF_ZERO);
201 if (m_field.check_field(nodals_positions)) {
202 CHKERR m_field.modify_finite_element_add_field_data(fe, nodals_positions);
203 }
204 CHKERR m_field.add_ents_to_finite_element_by_type(tris, MBTRI, fe);
206}
207
209 string problem) {
211 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(problem, ROW, &F);
212 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(problem, COL, &D);
214 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(problem, &A);
215
216 CHKERR KSPCreate(PETSC_COMM_WORLD, &kspSolver);
217 CHKERR KSPSetOperators(kspSolver, A, A);
218 CHKERR KSPSetFromOptions(kspSolver);
220}
221
223 string problem, string fe,
224 DirichletBC &bc,
225 Range &tris) {
227
228 CHKERR VecZeroEntries(F);
229 CHKERR MatZeroEntries(A);
230
233 CHKERR m_field.loop_finite_elements(problem, fe,
235 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
236 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
237 CHKERR VecAssemblyBegin(F);
238 CHKERR VecAssemblyEnd(F);
239 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
240 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
241
242 CHKERR KSPSolve(kspSolver, F, D);
243 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
244 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
245
246 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
247 problem, ROW, D, INSERT_VALUES, SCATTER_REVERSE);
248
249 bc.trisPtr = boost::shared_ptr<Range>(new Range(tris));
250 bc.mapZeroRows.clear();
251 bc.dofsIndices.clear();
252 bc.dofsValues.clear();
253
255}
256
258 string problem, string fe,
259 DirichletBC &bc) {
261 EntityHandle fe_meshset = m_field.get_finite_element_meshset("BC_FE");
262 Range bc_tris;
263 CHKERR m_field.get_moab().get_entities_by_type(fe_meshset, MBTRI, bc_tris);
264 return solveProblem(m_field, problem, fe, bc, bc_tris);
266}
267
270 CHKERR KSPDestroy(&kspSolver);
271 CHKERR MatDestroy(&A);
272 CHKERR VecDestroy(&F);
273 CHKERR VecDestroy(&D);
275}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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 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.
Definition: Exceptions.hpp:56
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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
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.
Structure used to enforce analytical boundary conditions.
DirichletBC(MoFEM::Interface &m_field, const std::string &field, Mat A, Vec X, Vec F)
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.
Definition: DirichletBC.hpp:57
std::map< DofIdx, FieldData > mapZeroRows
Definition: DirichletBC.hpp:70
std::vector< int > dofsIndices
Definition: DirichletBC.hpp:71
std::vector< double > dofsValues
Definition: DirichletBC.hpp:72
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.
Vec & snes_f
residual
Mat & snes_B
preconditioner of jacobian matrix
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23