OOFEM 3.0
Loading...
Searching...
No Matches
mixedgradientpressureneumann.C
Go to the documentation of this file.
1/*
2 *
3 * ##### ##### ###### ###### ### ###
4 * ## ## ## ## ## ## ## ### ##
5 * ## ## ## ## #### #### ## # ##
6 * ## ## ## ## ## ## ## ##
7 * ## ## ## ## ## ## ## ##
8 * ##### ##### ## ###### ## ##
9 *
10 *
11 * OOFEM : Object Oriented Finite Element Code
12 *
13 * Copyright (C) 1993 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "mixedgradientpressureneumann.h"
36#include "dofiditem.h"
37#include "dofmanager.h"
38#include "dof.h"
39#include "valuemodetype.h"
40#include "floatarray.h"
41#include "floatmatrix.h"
42#include "engngm.h"
43#include "node.h"
44#include "element.h"
45#include "integrationrule.h"
46#include "gaussintegrationrule.h"
47#include "gausspoint.h"
48#include "masterdof.h"
49#include "classfactory.h" // For sparse matrix creation.
50#include "sparsemtrxtype.h"
51#include "mathfem.h"
52#include "sparsemtrx.h"
53#include "sparselinsystemnm.h"
54#include "set.h"
55#include "dynamicinputrecord.h"
56#include "feinterpol.h"
57#include "unknownnumberingscheme.h"
58
59#ifdef _OPENMP
60#include <omp.h>
61#endif
62
63namespace oofem {
64REGISTER_BoundaryCondition(MixedGradientPressureNeumann);
65
66MixedGradientPressureNeumann :: MixedGradientPressureNeumann(int n, Domain *d) : MixedGradientPressureBC(n, d),
67 sigmaDev(std::make_unique<Node>(-1, d)) // Node number lacks meaning here.
68{
70 int nsd = d->giveNumberOfSpatialDimensions();
71 int components = nsd * nsd - 1;
72 for ( int i = 0; i < components; i++ ) {
73 int dofid = d->giveNextFreeDofID();
74 dev_id.followedBy(dofid);
75 // Just putting in X_i id-items since they don't matter.
76 sigmaDev->appendDof( new MasterDof( sigmaDev.get(), ( DofIDItem )dofid ) );
77 }
78}
79
80
81MixedGradientPressureNeumann :: ~MixedGradientPressureNeumann()
82{
83}
84
85
86int MixedGradientPressureNeumann :: giveNumberOfInternalDofManagers()
87{
88 return 1;
89}
90
91
92DofManager *MixedGradientPressureNeumann :: giveInternalDofManager(int i)
93{
94 return this->sigmaDev.get();
95}
96
97
98void MixedGradientPressureNeumann :: fromDeviatoricBase2D(FloatArray &cartesian, FloatArray &deviatoric)
99{
100 cartesian.resize(3);
101 cartesian.at(1) = deviatoric.at(1) / sqrt(2.0);
102 cartesian.at(2) = -deviatoric.at(1) / sqrt(2.0);
103 cartesian.at(3) = ( deviatoric.at(2) + deviatoric.at(3) ) * 0.5;
104}
105
106
107void MixedGradientPressureNeumann :: fromDeviatoricBase3D(FloatArray &cartesian, FloatArray &deviatoric)
108{
109 double s6 = sqrt(6.);
110 double s2 = sqrt(2.);
111 cartesian.resize(6);
112 cartesian.at(1) = 2.0 * deviatoric.at(1) / s6;
113 cartesian.at(2) = -deviatoric.at(1) / s6 + deviatoric.at(2) / s2;
114 cartesian.at(3) = -deviatoric.at(1) / s6 - deviatoric.at(2) / s2;
115 //
116 cartesian.at(4) = ( deviatoric.at(3) + deviatoric.at(6) ) * 0.5;
117 cartesian.at(5) = ( deviatoric.at(4) + deviatoric.at(7) ) * 0.5;
118 cartesian.at(6) = ( deviatoric.at(5) + deviatoric.at(8) ) * 0.5;
119}
120
121
122void MixedGradientPressureNeumann :: fromDeviatoricBase2D(FloatMatrix &cartesian, FloatMatrix &deviatoric)
123{
124 cartesian.resize(3, 3);
125 // E1 = [1/sqrt(2), 1/sqrt(2), 0,0]'; E2 = [0,0,1,0]'; E3 = [0,0,0,1]';
126 // C = E1*(E1'*e11 + E2'*e12 + E3'*e13) + E2*(E1'*e21 + E2'*e22 + E3'*e23) + E3*(E1'*e31 + E2'*e32 + E3'*e33);
127
128 cartesian.at(1, 1) = deviatoric.at(1, 1) * 0.5;
129 cartesian.at(2, 2) = deviatoric.at(1, 1) * 0.5;
130 cartesian.at(1, 2) = -deviatoric.at(1, 1) * 0.5;
131 cartesian.at(2, 1) = -deviatoric.at(1, 1) * 0.5;
132 //
133 cartesian.at(1, 3) = ( deviatoric.at(1, 2) + deviatoric.at(1, 3) ) / sqrt(8.0);
134 cartesian.at(2, 3) = -( deviatoric.at(1, 2) + deviatoric.at(1, 3) ) / sqrt(8.0);
135 cartesian.at(3, 1) = ( deviatoric.at(2, 1) + deviatoric.at(3, 1) ) / sqrt(8.0);
136 cartesian.at(3, 2) = -( deviatoric.at(2, 1) + deviatoric.at(3, 1) ) / sqrt(8.0);
137 //
138 cartesian.at(3, 3) = ( deviatoric.at(2, 2) + deviatoric.at(2, 3) +
139 deviatoric.at(3, 2) + deviatoric.at(3, 3) ) * 0.25;
140}
141
142
143void MixedGradientPressureNeumann :: fromDeviatoricBase3D(FloatMatrix &cartesian, FloatMatrix &deviatoric)
144{
145 cartesian.resize(6, 6);
146 /*
147 * syms
148 * e11 e12 e13 e14 e15 e16 e17 e18
149 * e21 e22 e23 e24 e25 e26 e27 e28
150 * e31 e32 e33 e34 e35 e36 e37 e38
151 * e41 e42 e43 e44 e45 e46 e47 e48
152 * e51 e52 e53 e54 e55 e56 e57 e58
153 * e61 e62 e63 e64 e65 e66 e67 e68
154 * e71 e72 e73 e74 e75 e76 e77 e78
155 * e81 e82 e83 e84 e85 e86 e87 e88;
156 * C = ...
157 * E1*(E1'*e11 + E2'*e12 + E3'*e13 + E4'*e14 + E5'*e15 + E6'*e16 + E7'*e17 + E8'*e18) + ...
158 * E2*(E1'*e21 + E2'*e22 + E3'*e23 + E4'*e24 + E5'*e25 + E6'*e26 + E7'*e27 + E8'*e28) + ...
159 * E3*(E1'*e31 + E2'*e32 + E3'*e33 + E4'*e34 + E5'*e35 + E6'*e36 + E7'*e37 + E8'*e38) + ...
160 * E4*(E1'*e41 + E2'*e42 + E3'*e43 + E4'*e44 + E5'*e45 + E6'*e46 + E7'*e47 + E8'*e48) + ...
161 * E5*(E1'*e51 + E2'*e52 + E3'*e53 + E4'*e54 + E5'*e55 + E6'*e56 + E7'*e57 + E8'*e58) + ...
162 * E6*(E1'*e61 + E2'*e62 + E3'*e63 + E4'*e64 + E5'*e65 + E6'*e66 + E7'*e67 + E8'*e68) + ...
163 * E7*(E1'*e71 + E2'*e72 + E3'*e73 + E4'*e74 + E5'*e75 + E6'*e76 + E7'*e77 + E8'*e78) + ...
164 * E8*(E1'*e81 + E2'*e82 + E3'*e83 + E4'*e84 + E5'*e85 + E6'*e86 + E7'*e87 + E8'*e88);
165 */
166
167 cartesian.at(1, 1) = deviatoric.at(1, 1) * 2.0 / 3.0;
168 cartesian.at(1, 2) = -deviatoric.at(1, 1) / 3.0 + deviatoric.at(1, 2) * sqrt(3.0);
169 cartesian.at(2, 1) = -deviatoric.at(1, 1) / 3.0 + deviatoric.at(1, 2) * sqrt(3.0);
170 cartesian.at(1, 3) = -deviatoric.at(1, 1) / 3.0 - deviatoric.at(1, 2) * sqrt(3.0);
171 cartesian.at(3, 1) = -deviatoric.at(1, 1) / 3.0 - deviatoric.at(1, 2) * sqrt(3.0);
172
173 cartesian.at(2, 2) = deviatoric.at(1, 1) / 6.0 + deviatoric.at(2, 2) / 2.0 - deviatoric.at(1, 2) / sqrt(12.0) - deviatoric.at(2, 1) / sqrt(12.0);
174 cartesian.at(2, 3) = deviatoric.at(1, 1) / 6.0 - deviatoric.at(2, 2) / 2.0 + deviatoric.at(1, 2) / sqrt(12.0) - deviatoric.at(2, 1) / sqrt(12.0);
175 cartesian.at(3, 3) = deviatoric.at(1, 1) / 6.0 + deviatoric.at(2, 2) / 2.0 - deviatoric.at(1, 2) / sqrt(12.0) + deviatoric.at(2, 1) / sqrt(12.0);
176 cartesian.at(3, 2) = deviatoric.at(1, 1) / 6.0 - deviatoric.at(2, 2) / 2.0 + deviatoric.at(1, 2) / sqrt(12.0) + deviatoric.at(2, 1) / sqrt(12.0);
177 // upper off diagonal part
178 cartesian.at(1, 4) = ( deviatoric.at(1, 3) + deviatoric.at(1, 6) ) / sqrt(6.0);
179 cartesian.at(1, 5) = ( deviatoric.at(1, 4) + deviatoric.at(1, 7) ) / sqrt(6.0);
180 cartesian.at(1, 6) = ( deviatoric.at(1, 5) + deviatoric.at(1, 8) ) / sqrt(6.0);
181 cartesian.at(2, 4) = deviatoric.at(2, 3) / sqrt(8.0) + deviatoric.at(2, 6) / sqrt(8.0) - deviatoric.at(1, 3) / sqrt(24.0) - deviatoric.at(1, 6) / sqrt(24.0);
182 cartesian.at(2, 5) = deviatoric.at(2, 4) / sqrt(8.0) + deviatoric.at(2, 7) / sqrt(8.0) - deviatoric.at(1, 4) / sqrt(24.0) - deviatoric.at(1, 7) / sqrt(24.0);
183 cartesian.at(2, 6) = deviatoric.at(2, 5) / sqrt(8.0) + deviatoric.at(2, 8) / sqrt(8.0) - deviatoric.at(1, 5) / sqrt(24.0) - deviatoric.at(1, 8) / sqrt(24.0);
184 cartesian.at(3, 4) = -deviatoric.at(2, 3) / sqrt(8.0) - deviatoric.at(2, 6) / sqrt(8.0) - deviatoric.at(1, 3) / sqrt(24.0) - deviatoric.at(1, 6) / sqrt(24.0);
185 cartesian.at(3, 5) = -deviatoric.at(2, 4) / sqrt(8.0) - deviatoric.at(2, 7) / sqrt(8.0) - deviatoric.at(1, 4) / sqrt(24.0) - deviatoric.at(1, 7) / sqrt(24.0);
186 cartesian.at(3, 6) = -deviatoric.at(2, 5) / sqrt(8.0) - deviatoric.at(2, 8) / sqrt(8.0) - deviatoric.at(1, 5) / sqrt(24.0) - deviatoric.at(1, 8) / sqrt(24.0);
187 // lower off diagonal part
188 cartesian.at(1, 4) = ( deviatoric.at(3, 1) + deviatoric.at(6, 1) ) / sqrt(6.0);
189 cartesian.at(1, 5) = ( deviatoric.at(4, 1) + deviatoric.at(7, 1) ) / sqrt(6.0);
190 cartesian.at(1, 6) = ( deviatoric.at(5, 1) + deviatoric.at(8, 1) ) / sqrt(6.0);
191 cartesian.at(2, 4) = deviatoric.at(3, 2) / sqrt(8.0) + deviatoric.at(6, 2) / sqrt(8.0) - deviatoric.at(3, 1) / sqrt(24.0) - deviatoric.at(6, 1) / sqrt(24.0);
192 cartesian.at(2, 5) = deviatoric.at(4, 2) / sqrt(8.0) + deviatoric.at(7, 2) / sqrt(8.0) - deviatoric.at(4, 1) / sqrt(24.0) - deviatoric.at(7, 1) / sqrt(24.0);
193 cartesian.at(2, 6) = deviatoric.at(5, 2) / sqrt(8.0) + deviatoric.at(8, 2) / sqrt(8.0) - deviatoric.at(5, 1) / sqrt(24.0) - deviatoric.at(8, 1) / sqrt(24.0);
194 cartesian.at(3, 4) = -deviatoric.at(3, 2) / sqrt(8.0) - deviatoric.at(6, 2) / sqrt(8.0) - deviatoric.at(3, 1) / sqrt(24.0) - deviatoric.at(6, 1) / sqrt(24.0);
195 cartesian.at(3, 5) = -deviatoric.at(4, 2) / sqrt(8.0) - deviatoric.at(7, 2) / sqrt(8.0) - deviatoric.at(4, 1) / sqrt(24.0) - deviatoric.at(7, 1) / sqrt(24.0);
196 cartesian.at(3, 6) = -deviatoric.at(5, 2) / sqrt(8.0) - deviatoric.at(8, 2) / sqrt(8.0) - deviatoric.at(5, 1) / sqrt(24.0) - deviatoric.at(8, 1) / sqrt(24.0);
197 //
198 cartesian.at(4, 4) = ( deviatoric.at(3, 3) + deviatoric.at(3, 6) + deviatoric.at(6, 3) + deviatoric.at(6, 6) ) * 0.25;
199 cartesian.at(4, 5) = ( deviatoric.at(3, 4) + deviatoric.at(3, 7) + deviatoric.at(6, 4) + deviatoric.at(6, 7) ) * 0.25;
200 cartesian.at(4, 6) = ( deviatoric.at(3, 5) + deviatoric.at(3, 8) + deviatoric.at(6, 5) + deviatoric.at(6, 8) ) * 0.25;
201 cartesian.at(5, 4) = ( deviatoric.at(4, 3) + deviatoric.at(4, 6) + deviatoric.at(7, 3) + deviatoric.at(7, 6) ) * 0.25;
202 cartesian.at(5, 5) = ( deviatoric.at(4, 4) + deviatoric.at(4, 7) + deviatoric.at(7, 4) + deviatoric.at(7, 7) ) * 0.25;
203 cartesian.at(5, 6) = ( deviatoric.at(4, 5) + deviatoric.at(4, 8) + deviatoric.at(7, 5) + deviatoric.at(7, 8) ) * 0.25;
204 cartesian.at(6, 4) = ( deviatoric.at(5, 3) + deviatoric.at(5, 6) + deviatoric.at(8, 3) + deviatoric.at(8, 6) ) * 0.25;
205 cartesian.at(6, 5) = ( deviatoric.at(5, 4) + deviatoric.at(5, 7) + deviatoric.at(8, 4) + deviatoric.at(8, 7) ) * 0.25;
206 cartesian.at(6, 6) = ( deviatoric.at(5, 5) + deviatoric.at(5, 8) + deviatoric.at(8, 5) + deviatoric.at(8, 8) ) * 0.25;
207}
208
209
210void MixedGradientPressureNeumann :: setPrescribedDeviatoricGradientFromVoigt(const FloatArray &t)
211{
212 // Converts the Voigt vector to a deviatoric base dyads (which don't assume symmetry)
213 int nsd = this->domain->giveNumberOfSpatialDimensions();
214 if ( nsd == 3 ) {
215 this->devGradient.resize(8);
216 this->devGradient.at(1) = ( 2.0 * t.at(1) - t.at(2) - t.at(3) ) / sqrt(6.);
217 this->devGradient.at(2) = ( t.at(2) - t.at(3) ) / sqrt(2.);
218 //
219 this->devGradient.at(3) = 0.5 * t.at(4);
220 this->devGradient.at(4) = 0.5 * t.at(5);
221 this->devGradient.at(5) = 0.5 * t.at(6);
222 //
223 this->devGradient.at(6) = 0.5 * t.at(4);
224 this->devGradient.at(7) = 0.5 * t.at(5);
225 this->devGradient.at(8) = 0.5 * t.at(6);
226
227 this->volGradient = t.at(1) + t.at(2) + t.at(3);
228 } else if ( nsd == 2 ) {
229 this->devGradient.resize(3);
230 this->devGradient.at(1) = ( t.at(1) - t.at(2) ) / sqrt(2.);
231 this->devGradient.at(2) = 0.5 * t.at(3);
232 this->devGradient.at(3) = 0.5 * t.at(3);
233
234 this->volGradient = t.at(1) + t.at(2);
235 } else {
236 this->devGradient.clear();
237 this->volGradient = t.at(1);
238 }
239}
240
241
242void MixedGradientPressureNeumann :: giveLocationArrays(std :: vector< IntArray > &rows, std :: vector< IntArray > &cols, CharType type,
243 const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
244{
245 if ( type != TangentStiffnessMatrix ) {
246 return;
247 }
248
249 IntArray loc_r, loc_c, sigma_loc_r, sigma_loc_c;
250
251 // Fetch the columns/rows for the stress contributions;
252 this->sigmaDev->giveLocationArray(dev_id, sigma_loc_r, r_s);
253 this->sigmaDev->giveLocationArray(dev_id, sigma_loc_c, c_s);
254
255 Set *set = this->giveDomain()->giveSet(this->set);
256 const IntArray &boundaries = set->giveBoundaryList();
257
258 rows.resize( boundaries.giveSize() );
259 cols.resize( boundaries.giveSize() );
260 int i = 0;
261 for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
262 Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
263 int boundary = boundaries.at(pos * 2);
264
265 const auto &bNodes = e->giveInterpolation()->boundaryGiveNodes(boundary, e->giveGeometryType());
266 e->giveBoundaryLocationArray(loc_r, bNodes, this->dofs, r_s);
267 e->giveBoundaryLocationArray(loc_c, bNodes, this->dofs, c_s);
268 // For most uses, loc_r == loc_c, and sigma_loc_r == sigma_loc_c.
269 rows [ i ] = loc_r;
270 cols [ i ] = sigma_loc_c;
271 i++;
272 // and the symmetric part (usually the transpose of above)
273 rows [ i ] = sigma_loc_r;
274 cols [ i ] = loc_c;
275 i++;
276 }
277}
278
279
280void MixedGradientPressureNeumann :: integrateVolTangent(FloatArray &answer, Element *e, int boundary)
281{
282 FloatArray normal, n;
283 FloatMatrix nMatrix;
284
285 FEInterpolation *interp = e->giveInterpolation(); // Geometry interpolation
286 // Some assumptions here, either velocity or displacement unknowns. Perhaps its better to just as for a certain equation id?
287 // Code will be written using v for velocity, but it could represent displacements.
288 FEInterpolation *interpUnknown = e->giveInterpolation(V_u);
289 if ( interpUnknown ) {
290 interpUnknown = e->giveInterpolation(D_u);
291 }
292
293 int nsd = e->giveDomain()->giveNumberOfSpatialDimensions();
294 // Order here should be the normal (which takes the first derivative) thus -1
295 int order = interp->giveInterpolationOrder() - 1 + interpUnknown->giveInterpolationOrder();
296 std :: unique_ptr< IntegrationRule > ir( interp->giveBoundaryIntegrationRule(order, boundary, e->giveGeometryType()) );
297
298 answer.clear();
299 for ( GaussPoint *gp: *ir ) {
300 const FloatArray &lcoords = gp->giveNaturalCoordinates();
301 FEIElementGeometryWrapper cellgeo(e);
302
303 // Evaluate the normal;
304 double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
305 // Evaluate the velocity/displacement coefficients
306 interpUnknown->boundaryEvalN(n, boundary, lcoords, cellgeo);
307 nMatrix.beNMatrixOf(n, nsd);
308
309 answer.plusProduct( nMatrix, normal, detJ * gp->giveWeight() );
310 }
311}
312
313
314void MixedGradientPressureNeumann :: integrateDevTangent(FloatMatrix &answer, Element *e, int boundary)
315{
316 FloatArray normal, n;
317 FloatMatrix nMatrix, E_n;
318 FloatMatrix contrib;
319
320 FEInterpolation *interp = e->giveInterpolation(); // Geometry interpolation
321 // Some assumptions here, either velocity or displacement unknowns. Perhaps its better to just as for a certain equation id?
322 // Code will be written using v for velocity, but it could represent displacements.
323 FEInterpolation *interpUnknown = e->giveInterpolation(V_u);
324 if ( interpUnknown ) {
325 interpUnknown = e->giveInterpolation(D_u);
326 }
327
328 int nsd = e->giveDomain()->giveNumberOfSpatialDimensions();
329 // Order here should be the normal (which takes the first derivative) thus -1
330 int order = interp->giveInterpolationOrder() - 1 + interpUnknown->giveInterpolationOrder();
331 std :: unique_ptr< IntegrationRule > ir( interp->giveBoundaryIntegrationRule(order, boundary, e->giveGeometryType()) );
332
333 answer.clear();
334 for ( auto &gp: *ir ) {
335 const FloatArray &lcoords = gp->giveNaturalCoordinates();
336 FEIElementGeometryWrapper cellgeo(e);
337
338 // Evaluate the normal;
339 double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
340 // Evaluate the velocity/displacement coefficients
341 interpUnknown->boundaryEvalN(n, boundary, lcoords, cellgeo);
342 nMatrix.beNMatrixOf(n, nsd);
343
344 // Formulating like this to avoid third order tensors, which is hard to express in linear algebra.
345 // v (x) n : E_i = v . ( E_i . n ) = v . E_n
346 if ( nsd == 3 ) {
347 E_n.resize(8, 3);
348 E_n.at(1, 1) = 2.0 * normal.at(1) / sqrt(6.0);
349 E_n.at(1, 2) = -normal.at(2) / sqrt(6.0);
350 E_n.at(1, 3) = -normal.at(3) / sqrt(6.0);
351
352 E_n.at(2, 1) = 0.;
353 E_n.at(2, 2) = normal.at(2) / sqrt(2.0);
354 E_n.at(2, 3) = -normal.at(3) / sqrt(2.0);
355
356 E_n.at(3, 1) = normal.at(2);
357 E_n.at(3, 2) = 0.;
358 E_n.at(3, 3) = 0.;
359
360 E_n.at(4, 1) = normal.at(3);
361 E_n.at(4, 2) = 0.;
362 E_n.at(4, 3) = 0.;
363
364 E_n.at(5, 1) = 0.;
365 E_n.at(5, 2) = normal.at(3);
366 E_n.at(5, 3) = 0.;
367
368 E_n.at(6, 1) = 0.;
369 E_n.at(6, 2) = normal.at(1);
370 E_n.at(6, 3) = 0.;
371
372 E_n.at(7, 1) = 0.;
373 E_n.at(7, 2) = 0.;
374 E_n.at(7, 3) = normal.at(1);
375
376 E_n.at(8, 1) = 0.;
377 E_n.at(8, 2) = 0.;
378 E_n.at(8, 3) = normal.at(2);
379 } else if ( nsd == 2 ) {
380 E_n.resize(3, 2);
381 E_n.at(1, 1) = normal.at(1) / sqrt(2.0);
382 E_n.at(1, 2) = -normal.at(2) / sqrt(2.0);
383
384 E_n.at(2, 1) = normal.at(2);
385 E_n.at(2, 2) = 0.;
386
387 E_n.at(3, 1) = 0.;
388 E_n.at(3, 2) = normal.at(1);
389 } else {
390 E_n.clear();
391 }
392
393 contrib.beProductOf(E_n, nMatrix);
394
395 answer.add(detJ * gp->giveWeight(), contrib);
396 }
397}
398
399
400void MixedGradientPressureNeumann :: assembleVector(FloatArray &answer, TimeStep *tStep,
401 CharType type, ValueModeType mode,
402 const UnknownNumberingScheme &s,
403 FloatArray *eNorms,
404 void* lock)
405{
406 Set *set = this->giveDomain()->giveSet(this->set);
407 const IntArray &boundaries = set->giveBoundaryList();
408
409 IntArray loc, sigma_loc; // For the velocities and stress respectively
410 IntArray masterDofIDs;
411 this->sigmaDev->giveLocationArray(dev_id, sigma_loc, s);
412
413 if ( type == ExternalForcesVector ) {
414 // The external forces have two contributions. On the additional equations for sigmaDev, the load is simple the deviatoric gradient.
415 double rve_size = this->domainSize();
416 FloatArray devLoad;
417 devLoad.beScaled(-rve_size, this->devGradient);
418#ifdef _OPENMP
419 if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
420#endif
421 answer.assemble(devLoad, sigma_loc);
422#ifdef _OPENMP
423 if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
424#endif
425
426 // The second contribution is on the momentumbalance equation; - int delta_v . n dA * p
427 FloatArray fe;
428 for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
429 Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
430 int boundary = boundaries.at(pos * 2);
431
432 const auto &bNodes = e->giveInterpolation()->boundaryGiveNodes(boundary, e->giveGeometryType());
433 e->giveBoundaryLocationArray(loc, bNodes, this->dofs, s, & masterDofIDs);
434 this->integrateVolTangent(fe, e, boundary);
435 fe.times(-this->pressure);
436#ifdef _OPENMP
437 if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
438#endif
439 answer.assemble(fe, loc);
440 if ( eNorms ) {
441 eNorms->assembleSquared(fe, masterDofIDs);
442 }
443#ifdef _OPENMP
444 if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
445#endif
446 }
447 } else if ( type == InternalForcesVector ) {
448 FloatMatrix Ke;
449 FloatArray fe_v, fe_s;
450 FloatArray s_dev, e_v;
451
452 // Fetch the current values of the stress;
453 this->sigmaDev->giveUnknownVector(s_dev, dev_id, mode, tStep);
454
455 // Assemble: int delta_v (x) n dA : E_i s_i
456 // int v (x) n dA : E_i delta_s_i
457 for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
458 Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
459 int boundary = boundaries.at(pos * 2);
460
461 // Fetch the element information;
462 const auto &bNodes = e->giveInterpolation()->boundaryGiveNodes(boundary, e->giveGeometryType());
463 e->giveBoundaryLocationArray(loc, bNodes, this->dofs, s, & masterDofIDs);
464 e->computeBoundaryVectorOf(bNodes, this->dofs, mode, tStep, e_v);
465 this->integrateDevTangent(Ke, e, boundary);
466
467 // We just use the tangent, less duplicated code (the addition of sigmaDev is linear).
468 fe_v.beProductOf(Ke, e_v);
469 fe_s.beTProductOf(Ke, s_dev);
470 // Note: The terms appear negative in the equations:
471 fe_v.negated();
472 fe_s.negated();
473
474#ifdef _OPENMP
475 if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
476#endif
477 answer.assemble(fe_s, loc); // Contributions to delta_v equations
478 answer.assemble(fe_v, sigma_loc); // Contribution to delta_s_i equations
479 if ( eNorms ) {
480 eNorms->assembleSquared(fe_s, masterDofIDs);
481 eNorms->assembleSquared(fe_v, dev_id);
482 }
483#ifdef _OPENMP
484 if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
485#endif
486 }
487 }
488}
489
490
491void MixedGradientPressureNeumann :: assemble(SparseMtrx &answer, TimeStep *tStep,
492 CharType type, const UnknownNumberingScheme &r_s,
493 const UnknownNumberingScheme &c_s, double scale,
494 void* lock)
495{
496 if ( type == TangentStiffnessMatrix || type == SecantStiffnessMatrix || type == ElasticStiffnessMatrix ) {
497 FloatMatrix Ke, KeT;
498 IntArray loc_r, loc_c, sigma_loc_r, sigma_loc_c;
499 Set *set = this->giveDomain()->giveSet(this->set);
500 const IntArray &boundaries = set->giveBoundaryList();
501
502 // Fetch the columns/rows for the stress contributions;
503 this->sigmaDev->giveLocationArray(dev_id, sigma_loc_r, r_s);
504 this->sigmaDev->giveLocationArray(dev_id, sigma_loc_c, c_s);
505
506 for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
507 Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
508 int boundary = boundaries.at(pos * 2);
509
510 // Fetch the element information;
511 const auto &bNodes = e->giveInterpolation()->boundaryGiveNodes(boundary, e->giveGeometryType());
512 e->giveBoundaryLocationArray(loc_r, bNodes, this->dofs, r_s);
513 e->giveBoundaryLocationArray(loc_c, bNodes, this->dofs, c_s);
514 this->integrateDevTangent(Ke, e, boundary);
515 Ke.negated();
516 Ke.times(scale);
517 KeT.beTranspositionOf(Ke);
518#ifdef _OPENMP
519 if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
520#endif
521 answer.assemble(sigma_loc_r, loc_c, Ke); // Contribution to delta_s_i equations
522 answer.assemble(loc_r, sigma_loc_c, KeT); // Contributions to delta_v equations
523#ifdef _OPENMP
524 if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
525#endif
526 }
527 }
528}
529
530
531void MixedGradientPressureNeumann :: computeFields(FloatArray &sigmaDev, double &vol, TimeStep *tStep)
532{
533 int nsd = this->giveDomain()->giveNumberOfSpatialDimensions();
534 FloatArray sigmaDevBase;
535 Set *set = this->giveDomain()->giveSet(this->set);
536 const IntArray &boundaries = set->giveBoundaryList();
537
538 // Fetch the current values of the stress in deviatoric base;
539 sigmaDevBase.resize( this->sigmaDev->giveNumberOfDofs() );
540 for ( int i = 1; i <= this->sigmaDev->giveNumberOfDofs(); i++ ) {
541 sigmaDevBase.at(i) = this->sigmaDev->giveDofWithID(dev_id.at(i))->giveUnknown(VM_Total, tStep);
542 }
543 // Convert it back from deviatoric base:
544 if ( nsd == 3 ) {
545 this->fromDeviatoricBase3D(sigmaDev, sigmaDevBase);
546 } else if ( nsd == 2 ) {
547 this->fromDeviatoricBase2D(sigmaDev, sigmaDevBase);
548 } else {
549 sigmaDev.clear();
550 }
551
552 // Postprocessing; vol = int v . n dA
553 FloatArray unknowns, fe;
554 vol = 0.;
555 for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
556 Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
557 int boundary = boundaries.at(pos * 2);
558
559 const auto &bNodes = e->giveInterpolation()->boundaryGiveNodes(boundary, e->giveGeometryType());
560 e->computeBoundaryVectorOf(bNodes, this->dofs, VM_Total, tStep, unknowns);
561 this->integrateVolTangent(fe, e, boundary);
562 vol += fe.dotProduct(unknowns);
563 }
564 double rve_size = this->domainSize();
565 vol /= rve_size;
566 vol -= volGradient; // This is needed for consistency; We return the volumetric "residual" if a gradient with volumetric contribution is set.
567}
568
569
570void MixedGradientPressureNeumann :: computeTangents(FloatMatrix &Ed, FloatArray &Ep, FloatArray &Cd, double &Cp, TimeStep *tStep)
571{
572 // Fetch some information from the engineering model
573 EngngModel *rve = this->giveDomain()->giveEngngModel();
575 std :: unique_ptr< SparseLinearSystemNM > solver(
576 classFactory.createSparseLinSolver( ST_Petsc, this->domain, this->domain->giveEngngModel() ) ); // = rve->giveLinearSolver();
577 SparseMtrxType stype = solver->giveRecommendedMatrix(true);
578 EModelDefaultEquationNumbering fnum;
579 Set *set = this->giveDomain()->giveSet(this->set);
580 const IntArray &boundaries = set->giveBoundaryList();
581 double rve_size = this->domainSize();
582
583 // Set up and assemble tangent FE-matrix which will make up the sensitivity analysis for the macroscopic material tangent.
584 std :: unique_ptr< SparseMtrx > Kff( classFactory.createSparseMtrx(stype) );
585 if ( !Kff ) {
586 OOFEM_ERROR("Couldn't create sparse matrix of type %d\n", stype);
587 }
588 Kff->buildInternalStructure(rve, this->domain->giveNumber(), fnum);
589 rve->assemble(*Kff, tStep, TangentAssembler(TangentStiffness), fnum, this->domain);
590
591 // Setup up indices and locations
592 int neq = Kff->giveNumberOfRows();
593
594 // Indices and such of internal dofs
595 int ndev = this->sigmaDev->giveNumberOfDofs();
596
597 // Matrices and arrays for sensitivities
598 FloatMatrix ddev_pert(neq, ndev); // In fact, npeq should most likely equal ndev
599 FloatArray p_pert(neq); // RHS for d_dev [d_dev11, d_dev22, d_dev12] in 2D
600
601 FloatMatrix s_d(neq, ndev); // Sensitivity fields for d_dev
602 FloatArray s_p(neq); // Sensitivity fields for p
603
604 // Unit pertubations for d_dev
605 ddev_pert.zero();
606 for ( int i = 1; i <= ndev; ++i ) {
607 int eqn = this->sigmaDev->giveDofWithID(dev_id.at(i))->giveEquationNumber(fnum);
608 ddev_pert.at(eqn, i) = -1.0 * rve_size;
609 }
610
611 // Unit pertubation for d_p
612 p_pert.zero();
613 FloatArray fe;
614 IntArray loc;
615 for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
616 Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
617 int boundary = boundaries.at(pos * 2);
618
619 const auto &bNodes = e->giveInterpolation()->boundaryGiveNodes(boundary, e->giveGeometryType());
620 e->giveBoundaryLocationArray(loc, bNodes, this->dofs, fnum);
621 this->integrateVolTangent(fe, e, boundary);
622 fe.times(-1.0); // here d_p = 1.0
623 p_pert.assemble(fe, loc);
624 }
625
626 // Solve all sensitivities
627 solver->solve(*Kff, ddev_pert, s_d);
628 solver->solve(*Kff, p_pert, s_p);
629
630 // Extract the stress response from the solutions
631 FloatArray sigma_p(ndev);
632 FloatMatrix sigma_d(ndev, ndev);
633 for ( int i = 1; i <= ndev; ++i ) {
634 int eqn = this->sigmaDev->giveDofWithID(dev_id.at(i))->giveEquationNumber(fnum);
635 sigma_p.at(i) = s_p.at(eqn);
636 for ( int j = 1; j <= ndev; ++j ) {
637 sigma_d.at(i, j) = s_d.at(eqn, j);
638 }
639 }
640
641 // Post-process the volumetric rate of deformations in the sensitivity fields;
642 FloatArray e_d(ndev);
643 e_d.zero();
644 double e_p = 0.0;
645 for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
646 Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
647 int boundary = boundaries.at(pos * 2);
648
649 this->integrateVolTangent(fe, e, boundary);
650
651 const auto &bNodes = e->giveInterpolation()->boundaryGiveNodes(boundary, e->giveGeometryType());
652 e->giveBoundaryLocationArray(loc, bNodes, this->dofs, fnum);
653
654 // Using "loc" to pick out the relevant contributions. This won't work at all if there are local coordinate systems in these nodes
655 // or slave nodes etc. The goal is to compute the velocity from the sensitivity field, but we need to avoid going through the actual
656 // engineering model. If this ever becomes an issue it needs to perform the same steps as Element::giveUnknownVector does.
657 for ( int i = 1; i <= fe.giveSize(); ++i ) {
658 if ( loc.at(i) > 0 ) {
659 e_p += fe.at(i) * s_p.at( loc.at(i) );
660 for ( int j = 1; j <= ndev; ++j ) {
661 e_d.at(j) += fe.at(i) * s_d.at(loc.at(i), j);
662 }
663 }
664 }
665 }
666 e_p /= rve_size;
667 e_d.times(1. / rve_size);
668
669 // Now we need to express the tangents in the normal cartesian coordinate system (as opposed to the deviatoric base we use during computations
670 Cp = - e_p; // Scalar components are of course the same
671 double nsd = this->giveDomain()->giveNumberOfSpatialDimensions();
672 if ( nsd == 3 ) {
673 this->fromDeviatoricBase3D(Cd, e_d);
674 this->fromDeviatoricBase3D(Ep, sigma_p);
675 this->fromDeviatoricBase3D(Ed, sigma_d);
676 } else if ( nsd == 2 ) {
677 this->fromDeviatoricBase2D(Cd, e_d);
678 this->fromDeviatoricBase2D(Ep, sigma_p);
679 this->fromDeviatoricBase2D(Ed, sigma_d);
680 } else { // For 1D case, there simply are no deviatoric components!
681 Cd.clear();
682 Ep.clear();
683 Ed.clear();
684 }
685}
686
687
688void MixedGradientPressureNeumann :: initializeFrom(InputRecord &ir)
689{
690 MixedGradientPressureBC :: initializeFrom(ir);
691}
692
693
694void MixedGradientPressureNeumann :: giveInputRecord(DynamicInputRecord &input)
695{
696 MixedGradientPressureBC :: giveInputRecord(input);
697 input.setField(this->pressure, _IFT_MixedGradientPressure_pressure);
698 OOFEM_ERROR("Not supported yet");
699 //FloatArray devGradientVoigt;
700 //input.setField(devGradientVoigt, _IFT_MixedGradientPressure_devGradient);
701}
702
703
704
705void MixedGradientPressureNeumann :: scale(double s)
706{
707 devGradient.times(s);
708 pressure *= s;
709}
710} // end namespace oofem
REGISTER_BoundaryCondition
#define REGISTER_BoundaryCondition(class)
Definition classfactory.h:133
oofem::Domain::giveNextFreeDofID
int giveNextFreeDofID(int increment=1)
Definition domain.C:1519
oofem::Domain::giveNumberOfSpatialDimensions
int giveNumberOfSpatialDimensions()
Returns number of spatial dimensions.
Definition domain.C:1137
oofem::DynamicInputRecord::setField
void setField(int item, InputFieldType id)
Definition dynamicinputrecord.C:302
oofem::Element::giveBoundaryLocationArray
virtual void giveBoundaryLocationArray(IntArray &locationArray, const IntArray &bNodes, const UnknownNumberingScheme &s, IntArray *dofIds=NULL)
Definition element.C:485
oofem::Element::giveInterpolation
virtual FEInterpolation * giveInterpolation() const
Definition element.h:648
oofem::Element::computeBoundaryVectorOf
void computeBoundaryVectorOf(const IntArray &bNodes, const IntArray &dofIDMask, ValueModeType u, TimeStep *tStep, FloatArray &answer, bool padding=false)
Definition element.C:163
oofem::Element::giveGeometryType
virtual Element_Geometry_Type giveGeometryType() const =0
oofem::EngngModel::assemble
virtual void assemble(SparseMtrx &answer, TimeStep *tStep, const MatrixAssembler &ma, const UnknownNumberingScheme &s, Domain *domain)
Definition engngm.C:889
oofem::FEInterpolation::giveBoundaryIntegrationRule
virtual std::unique_ptr< IntegrationRule > giveBoundaryIntegrationRule(int order, int boundary, const Element_Geometry_Type) const
Definition feinterpol.C:101
oofem::FEInterpolation::boundaryGiveNodes
virtual IntArray boundaryGiveNodes(int boundary, const Element_Geometry_Type) const =0
oofem::FEInterpolation::giveInterpolationOrder
int giveInterpolationOrder() const
Definition feinterpol.h:199
oofem::FEInterpolation::boundaryEvalNormal
virtual double boundaryEvalNormal(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const =0
oofem::FEInterpolation::boundaryEvalN
virtual void boundaryEvalN(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const =0
oofem::FEMComponent::giveDomain
Domain * giveDomain() const
Definition femcmpnn.h:97
oofem::FEMComponent::domain
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition femcmpnn.h:79
oofem::FloatArray::assemble
void assemble(const FloatArray &fe, const IntArray &loc)
Definition floatarray.C:616
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::plusProduct
void plusProduct(const FloatMatrix &b, const FloatArray &s, double dV)
Definition floatarray.C:288
oofem::FloatArray::dotProduct
double dotProduct(const FloatArray &x) const
Definition floatarray.C:524
oofem::FloatArray::giveSize
Index giveSize() const
Returns the size of receiver.
Definition floatarray.h:261
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition floatarray.C:683
oofem::FloatArray::beProductOf
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Definition floatarray.C:689
oofem::FloatArray::assembleSquared
void assembleSquared(const FloatArray &fe, const IntArray &loc)
Definition floatarray.C:637
oofem::FloatArray::beScaled
void beScaled(double s, const FloatArray &b)
Definition floatarray.C:208
oofem::FloatArray::beTProductOf
void beTProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Definition floatarray.C:721
oofem::FloatArray::times
void times(double s)
Definition floatarray.C:834
oofem::FloatMatrix::times
void times(double f)
Definition floatmatrix.C:1155
oofem::FloatMatrix::add
void add(const FloatMatrix &a)
Definition floatmatrix.C:962
oofem::FloatMatrix::resize
void resize(Index rows, Index cols)
Definition floatmatrix.C:79
oofem::FloatMatrix::clear
*Sets size of receiver to be an empty matrix It will have zero rows and zero columns size void clear()
Definition floatmatrix.h:540
oofem::FloatMatrix::beProductOf
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Definition floatmatrix.C:365
oofem::FloatMatrix::beNMatrixOf
void beNMatrixOf(const FloatArray &n, int nsd)
Definition floatmatrix.C:1633
oofem::FloatMatrix::beTranspositionOf
void beTranspositionOf(const FloatMatrix &src)
Definition floatmatrix.C:351
oofem::FloatMatrix::zero
void zero()
Zeroes all coefficient of receiver.
Definition floatmatrix.C:1064
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::GeneralBoundaryCondition::set
int set
Set number for boundary condition to be applied to.
Definition generalboundarycondition.h:96
oofem::GeneralBoundaryCondition::dofs
IntArray dofs
Dofs that b.c. is applied to (relevant for Dirichlet type b.c.s).
Definition generalboundarycondition.h:89
oofem::IntArray::at
int & at(std::size_t i)
Definition intarray.h:104
oofem::IntArray::giveSize
int giveSize() const
Definition intarray.h:211
oofem::MixedGradientPressureNeumann::fromDeviatoricBase2D
void fromDeviatoricBase2D(FloatArray &cartesian, FloatArray &deviatoric)
Converts from deviatoric to (normal) cartesian base (arrays are second order 2D tensors in Voigt nota...
Definition mixedgradientpressureneumann.C:98
oofem::MixedGradientPressureNeumann::sigmaDev
std ::unique_ptr< Node > sigmaDev
DOF-manager containing the unknown deviatoric stress.
Definition mixedgradientpressureneumann.h:106
oofem::MixedGradientPressureNeumann::devGradient
FloatArray devGradient
Prescribed gradient in Voigt form.
Definition mixedgradientpressureneumann.h:95
oofem::MixedGradientPressureNeumann::integrateVolTangent
void integrateVolTangent(FloatArray &answer, Element *e, int boundary)
Helper function that integrates the volumetric tangent contribution from a single element boundary.
Definition mixedgradientpressureneumann.C:280
oofem::MixedGradientPressureNeumann::dev_id
IntArray dev_id
Dof IDs for the lagrange multipliers in sigmaDev.
Definition mixedgradientpressureneumann.h:108
oofem::MixedGradientPressureNeumann::fromDeviatoricBase3D
void fromDeviatoricBase3D(FloatArray &cartesian, FloatArray &deviatoric)
Converts from deviatoric to (normal) cartesian base (arrays are second order 3D tensors in Voigt nota...
Definition mixedgradientpressureneumann.C:107
oofem::MixedGradientPressureNeumann::integrateDevTangent
void integrateDevTangent(FloatMatrix &answer, Element *e, int boundary)
Helper function that integrates the deviatoric tangent contribution from a single element boundary.
Definition mixedgradientpressureneumann.C:314
oofem::SparseMtrx::assemble
virtual int assemble(const IntArray &loc, const FloatMatrix &mat)=0
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
_IFT_MixedGradientPressure_pressure
#define _IFT_MixedGradientPressure_pressure
Definition mixedgradientpressurebc.h:44
oofem::classFactory
ClassFactory & classFactory
Definition classfactory.C:92
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011