v0.13.1
BCs_RVELagrange_Trac.hpp
Go to the documentation of this file.
1/* Copyright (C) 2014, Zahur Ullah (Zahur.Ullah AT glasgow.ac.uk)
2 * --------------------------------------------------------------
3 */
4
5/* This file is part of MoFEM.
6 * MoFEM is free software: you can redistribute it and/or modify it under
7 * the terms of the GNU Lesser General Public License as published by the
8 * Free Software Foundation, either version 3 of the License, or (at your
9 * option) any later version.
10 *
11 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 * License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18
19#ifndef __BCS_RVELAGRANGE_TRAC_HPP
20#define __BCS_RVELAGRANGE_TRAC_HPP
21
23
25
26 struct CommonData {
28 };
30
32
33 PetscErrorCode shapeMat(
34 int rank,unsigned int gg,DataForcesAndSourcesCore::EntData &col_data,MatrixDouble &n
35 ) {
36 PetscFunctionBegin;
37
38 int shape_size = col_data.getN().size2();
39
40 n.resize(rank,shape_size*rank);
41 n.clear();
42
43 int kk = 0;
44 for(int ii=0; ii<shape_size; ii++){
45 double val = col_data.getN()(gg,ii);
46 for(int jj=0; jj<rank; jj++){
47 n(jj,kk) = val;
48 kk++;
49 }
50 }
51
52 PetscFunctionReturn(0);
53 }
54
55
57
58 PetscErrorCode hMat(
59 VectorDouble face_normal,
60 int rank,
61 MatrixDouble &n_mat,
62 MatrixDouble &h_mat
63 ) {
64 PetscFunctionBegin;
65
66 switch(rank) {
67 //Mechanical problem
68 case 3: {
69 H_mat_1Node.resize(6,3);
70 H_mat_1Node.clear(); //for one node
71 if(face_normal[0]>0) { //+X face of the RVE
72 H_mat_1Node(0,0) = 1.0;
73 H_mat_1Node(3,1) = 1.0;
74 H_mat_1Node(5,2) = 1.0;
75 } else if(face_normal[0]<0) { //-X face of the RVE
76 H_mat_1Node(0,0) = -1.0;
77 H_mat_1Node(3,1) = -1.0;
78 H_mat_1Node(5,2) = -1.0;
79 } else if(face_normal[1]>0) { //+Y face of the RVE
80 H_mat_1Node(1,1) = 1.0;
81 H_mat_1Node(3,0) = 1.0;
82 H_mat_1Node(4,2) = 1.0;
83 } else if(face_normal[1]<0) { //-Y face of the RVE
84 H_mat_1Node(1,1) = -1.0;
85 H_mat_1Node(3,0) = -1.0;
86 H_mat_1Node(4,2) = -1.0;
87 } else if(face_normal[2]>0) { //+Z face of the RVE
88 H_mat_1Node(2,2) = 1.0;
89 H_mat_1Node(4,1) = 1.0;
90 H_mat_1Node(5,0) = 1.0;
91 } else if(face_normal[2]<0) { //-Z face of the RVE
92 H_mat_1Node(2,2) = -1.0;
93 H_mat_1Node(4,1) = -1.0;
94 H_mat_1Node(5,0) = -1.0;
95 } else {
96 SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Can not be ?!");
97 }
98
99 int num_col = n_mat.size2();
100 h_mat.resize(6,num_col);
101
102 int cc1 = 0;
103 for(int bb = 0; bb<num_col/3; bb++) {
104 //blocks of 6x3
105 for(int rr = 0; rr<6; rr++) {
106 for(int cc = 0; cc<3; cc++) {
107 h_mat(rr,(cc+cc1)) = H_mat_1Node(rr,cc);
108 }
109 }
110 cc1+=3;
111 }
112
113 }
114 break;
115 case 1: //Moisture transport or thermal problem
116 {
117 H_mat_1Node.resize(3,1);
118 H_mat_1Node.clear();
119 if(face_normal[0]>0) { //+X face of the RVE
120 H_mat_1Node(0,0) = 1.0;
121 }
122 if(face_normal[0]<0) { //-X face of the RVE
123 H_mat_1Node(0,0) = -1.0;
124 }
125 if(face_normal[1]>0) { //+Y face of the RVE
126 H_mat_1Node(1,0) = 1.0;
127 }
128 if(face_normal[1]<0) { //-Y face of the RVE
129 H_mat_1Node(1,0) = -1.0;
130 }
131 if(face_normal[2]>0) { //+Z face of the RVE
132 H_mat_1Node(2,0) = 1.0;
133 }
134 if(face_normal[2]<0) { //-Z face of the RVE
135 H_mat_1Node(2,0) = -1.0;
136 }
137
138 int num_col = n_mat.size2();
139 h_mat.resize(3,num_col);
140
141 int cc1 = 0;
142 for(int bb = 0; bb<num_col; bb++) {
143 //blocks of 3x1
144 for(int rr = 0; rr<3; rr++) {
145 for(int cc = 0; cc<1; cc++) {
146 h_mat(rr,(cc+cc1))=H_mat_1Node(rr,cc);
147 }
148 }
149 cc1+=1;
150 }
151
152 }
153 break;
154 default:
155 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
156 }
157 PetscFunctionReturn(0);
158 }
159 };
160
162
163 /// \biref operator to calculate the LHS for the RVE bounary conditions
165
169 Mat Aij;
170
171 OpRVEBCsLhs(const string field_name, const string lagrang_field_name,Mat aij,
172 RVEBC_Data &data, CommonFunctions &common_functions, bool ho_geometry = false
173 ):
175 dAta(data),commonFunctions(common_functions),hoGeometry(ho_geometry),Aij(aij) {
176 //This will make sure to loop over all entities
177 // (e.g. for order=2 it will make doWork to loop 16 time)
178 sYmm = false;
179 }
180
182 ublas::vector<MatrixDouble > N_mat_row;
184
185 PetscErrorCode doWork(
186 int row_side,int col_side,
187 EntityType row_type,EntityType col_type,
190 ) {
191 PetscFunctionBegin;
192
193 try {
194 if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
195 if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
196
197 int nb_row = row_data.getIndices().size();
198 int nb_col = col_data.getIndices().size();
199
200 auto weak_ptr_dof =
201 getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
202 col_data.getIndices()[0]);
203 const FENumeredDofEntity *dof_ptr;
204 if (auto ptr = weak_ptr_dof.lock())
205 dof_ptr = ptr.get();
206 else
207 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
208
209 int rank = dof_ptr->getNbOfCoeffs();
210
211 K.resize(nb_row,nb_col);
212 K.clear();
213
214 N_mat_row.resize(col_data.getN().size1());
215 for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
216 double area;
217 if(hoGeometry) {
218 area = norm_2(getNormalsAtGaussPts(gg))*0.5;
219 } else {
220 area = getArea();
221 }
222 double val = getGaussPts()(2,gg)*area;
223
224 if(col_type==MBVERTEX) {
225 ierr = commonFunctions.shapeMat(rank,gg,col_data,N_mat_row[gg]); CHKERRQ(ierr);
226 }
227 ierr = commonFunctions.shapeMat(rank,gg,col_data,N_mat_col); CHKERRQ(ierr);
228 ierr = commonFunctions.hMat(getNormalsAtGaussPts(gg),rank,N_mat_row(gg),H_mat); CHKERRQ(ierr);
229
230 if(gg==0){
231 NTN = prod(trans(N_mat_row(gg)),N_mat_col); //we don't need to define size NTN
232 K = val*prod(H_mat,NTN); //we don't need to defiend size of K
233 } else {
234 noalias(NTN) = prod(trans(N_mat_row(gg)),N_mat_col);
235 K += val*prod(H_mat, NTN);
236 }
237 }
238
239 // Matrix C
240 int nb_rows=row_data.getIndices().size();
241 int nb_cols=col_data.getIndices().size();
243 Aij,nb_rows,&row_data.getIndices()[0],nb_cols,&col_data.getIndices()[0],&K(0,0),ADD_VALUES
244 ); CHKERRQ(ierr);
245
246 // Matrix CT
247 transK = trans(K);
249 Aij,nb_cols,&col_data.getIndices()[0],nb_rows,&row_data.getIndices()[0],&transK(0,0),ADD_VALUES
250 ); CHKERRQ(ierr);
251
252 } catch (const std::exception& ex) {
253 ostringstream ss;
254 ss << "throw in method: " << ex.what() << endl;
255 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,ss.str().c_str());
256 }
257 PetscFunctionReturn(0);
258 }
259
260 };
261
262 /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
264
269
271 const string field_name, RVEBC_Data &data,CommonData &common_data,CommonFunctions &common_functions,bool ho_geometry = false
272 ):
274 dAta(data),
275 hoGeometry(ho_geometry),
276 commonData(common_data),
278 }
279
281
283 PetscFunctionBegin;
284
285 if(type!=MBVERTEX) PetscFunctionReturn(0);
286 if(data.getIndices().size()==0) PetscFunctionReturn(0);
287
288 auto weak_ptr_dof =
289 getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
290 data.getIndices()[0]);
291 const FENumeredDofEntity *dof_ptr;
292 if (auto ptr = weak_ptr_dof.lock())
293 dof_ptr = ptr.get();
294 else
295 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
296
297 int rank = dof_ptr->getNbOfCoeffs();
298
299 double x,y,z;
300 for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
301 double area;
302 if(hoGeometry) {
303 area = norm_2(getNormalsAtGaussPts(gg))*0.5;
304 } else {
305 area = getArea();
306 }
307 double val = getGaussPts()(2,gg)*area;
308
309 x = getCoordsAtGaussPts()(gg,0);
310 y = getCoordsAtGaussPts()(gg,1);
311 z = getCoordsAtGaussPts()(gg,2);
312
313 switch(rank) {
314 case 3: //mech problem
315 X_mat.resize(rank,6,false);
316 X_mat.clear();
317 X_mat(0,0) = 2.0*x; X_mat(0,3) = y; X_mat(0,5) = z;
318 X_mat(1,1) = 2.0*y; X_mat(1,3) = x; X_mat(1,4) = z;
319 X_mat(2,2) = 2.0*z; X_mat(2,4) = y; X_mat(2,5) = x;
320 X_mat=0.5*X_mat;
321 break;
322 case 1: //moisture transport or thermal problem
323 X_mat.resize(rank,3,false);
324 X_mat.clear();
325 X_mat(0,0) = x;
326 X_mat(0,1) = y;
327 X_mat(0,2) = z;
328 break;
329 default:
330 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
331 }
332 ierr = commonFunctions.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);
333 ierr = commonFunctions.hMat(getNormalsAtGaussPts(gg), rank, N_mat, H_mat); CHKERRQ(ierr);
334
335 if(gg==0) {
336 NTX = val*prod(trans(N_mat),X_mat);
337 commonData.D_mat = prod(H_mat, NTX);
338 } else {
339 noalias(NTX) = val*prod(trans(N_mat),X_mat);
340 commonData.D_mat += prod(H_mat, NTX);
341 }
342 }
343
344 PetscFunctionReturn(0);
345 }
346 };
347
348 /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
350
353 vector<Vec> &F;
354
355 OpRVEBCsRhs_Assemble(const string lagrang_field_name,vector<Vec> &f,RVEBC_Data &data,CommonData &common_data):
357 dAta(data),commonData(common_data),F(f) {
358 }
359
360
362
364 PetscFunctionBegin;
365
366 if(data.getIndices().size()==0) PetscFunctionReturn(0);
367
369 applied_strain.resize(D_mat.size2(),false);
370
371 for(int ii = 0;ii<F.size();ii++) {
372 applied_strain.clear();
373 applied_strain(ii) = 1.0;
374 f.resize(D_mat.size1(),false);
375 noalias(f) = prod(D_mat,applied_strain);
376 ierr = VecSetValues(F[ii],data.getIndices().size(),&data.getIndices()[0],&f[0],ADD_VALUES); CHKERRQ(ierr);
377 }
378
379 PetscFunctionReturn(0);
380 }
381 };
382
383 PetscErrorCode setRVEBCsOperators(
384 string field_name,string lagrang_field_name,string mesh_nodals_positions,Mat aij,vector<Vec> &f
385 ) {
386 PetscFunctionBegin;
387 bool ho_geometry = false;
388 if(mField.check_field(mesh_nodals_positions)) {
389 ho_geometry = true;
390 }
391 map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
392 for(;sit!=setOfRVEBC.end();sit++) {
393 //LHS
394 feRVEBCLhs.getOpPtrVector().push_back(
395 new OpRVEBCsLhs(field_name,lagrang_field_name,aij,sit->second,commonFunctions,ho_geometry)
396 );
397 //RHS
398 feRVEBCRhs.getOpPtrVector().push_back(
399 new OpRVEBCsRhs_Cal(field_name,sit->second,commonData,commonFunctions,ho_geometry)
400 );
401 feRVEBCRhs.getOpPtrVector().push_back(
402 new OpRVEBCsRhs_Assemble(lagrang_field_name,f,sit->second,commonData)
403 );
404 }
405 PetscFunctionReturn(0);
406 }
407
409
413
415 const string lagrang_field_name,Vec stress_homo,RVEBC_Data &data,CommonData &common_data
416 ):
418 dAta(data),commonData(common_data),stressHomo(stress_homo) {
419 }
420
421
423
425 PetscFunctionBegin;
426 if(data.getIndices().size()==0) PetscFunctionReturn(0);
427 if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
428 stressVector.resize(commonData.D_mat.size2());
429 //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
430 noalias(stressVector) = prod(trans(commonData.D_mat),-data.getFieldData());
431 const int indices_6[6] = {0, 1, 2, 3, 4, 5};
432 const int indices_3[3] = {0, 1, 2};
433 switch(stressVector.size()) {
434 case 6:
435 ierr = VecSetValues(stressHomo,6,indices_6,&stressVector[0],ADD_VALUES); CHKERRQ(ierr);
436 break;
437 case 3:
438 ierr = VecSetValues(stressHomo,3,indices_3,&stressVector[0],ADD_VALUES); CHKERRQ(ierr);
439 break;
440 default:
441 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
442 }
443 PetscFunctionReturn(0);
444 }
445
446 };
447
449 string field_name,string lagrang_field_name,string mesh_nodals_positions,Vec stress_homo
450 ) {
451 PetscFunctionBegin;
452 bool ho_geometry = false;
453 if(mField.check_field(mesh_nodals_positions)) {
454 ho_geometry = true;
455 }
456 map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
457 for(;sit!=setOfRVEBC.end();sit++) {
458 feRVEBCStress.getOpPtrVector().push_back(
460 );
461 feRVEBCStress.getOpPtrVector().push_back(
462 new OpRVEHomoStress_Assemble(lagrang_field_name,stress_homo,sit->second,commonData)
463 );
464 }
465 PetscFunctionReturn(0);
466 }
467
468};
469
470#endif
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
FTensor::Index< 'n', SPACE_DIM > n
virtual bool check_field(const std::string &name) const =0
check if field is in database
const FTensor::Tensor2< T, Dim, Dim > Vec
auto f
Definition: HenckyOps.hpp:5
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
constexpr auto field_name
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
CommonFunctions common_functions
MoFEM::Interface & mField
PetscErrorCode hMat(VectorDouble face_normal, int rank, MatrixDouble &n_mat, MatrixDouble &h_mat)
PetscErrorCode shapeMat(int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &n)
\biref operator to calculate the LHS for the RVE bounary conditions
ublas::vector< MatrixDouble > N_mat_row
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpRVEBCsLhs(const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data &data, CommonFunctions &common_functions, bool ho_geometry=false)
\biref operator to calculate the RHS of the constrain for the RVE boundary conditions
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsRhs_Assemble(const string lagrang_field_name, vector< Vec > &f, RVEBC_Data &data, CommonData &common_data)
\biref operator to calculate the RHS of the constrain for the RVE boundary conditions
OpRVEBCsRhs_Cal(const string field_name, RVEBC_Data &data, CommonData &common_data, CommonFunctions &common_functions, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEHomoStress_Assemble(const string lagrang_field_name, Vec stress_homo, RVEBC_Data &data, CommonData &common_data)
BCs_RVELagrange_Trac(MoFEM::Interface &m_field)
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Deprecated interface functions.