OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
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 - 2013 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 namespace oofem {
60 REGISTER_BoundaryCondition(MixedGradientPressureNeumann);
61 
62 MixedGradientPressureNeumann :: MixedGradientPressureNeumann(int n, Domain *d) : MixedGradientPressureBC(n, d),
63  sigmaDev(new Node(-1, d)) // Node number lacks meaning here.
64 {
66  int nsd = d->giveNumberOfSpatialDimensions();
67  int components = nsd * nsd - 1;
68  for ( int i = 0; i < components; i++ ) {
69  int dofid = d->giveNextFreeDofID();
70  dev_id.followedBy(dofid);
71  // Just putting in X_i id-items since they don't matter.
72  sigmaDev->appendDof( new MasterDof( sigmaDev.get(), ( DofIDItem )dofid ) );
73  }
74 }
75 
76 
77 MixedGradientPressureNeumann :: ~MixedGradientPressureNeumann()
78 {
79 }
80 
81 
82 int MixedGradientPressureNeumann :: giveNumberOfInternalDofManagers()
83 {
84  return 1;
85 }
86 
87 
88 DofManager *MixedGradientPressureNeumann :: giveInternalDofManager(int i)
89 {
90  return this->sigmaDev.get();
91 }
92 
93 
94 void MixedGradientPressureNeumann :: fromDeviatoricBase2D(FloatArray &cartesian, FloatArray &deviatoric)
95 {
96  cartesian.resize(3);
97  cartesian.at(1) = deviatoric.at(1) / sqrt(2.0);
98  cartesian.at(2) = -deviatoric.at(1) / sqrt(2.0);
99  cartesian.at(3) = ( deviatoric.at(2) + deviatoric.at(3) ) * 0.5;
100 }
101 
102 
103 void MixedGradientPressureNeumann :: fromDeviatoricBase3D(FloatArray &cartesian, FloatArray &deviatoric)
104 {
105  double s6 = sqrt(6.);
106  double s2 = sqrt(2.);
107  cartesian.resize(6);
108  cartesian.at(1) = 2.0 * deviatoric.at(1) / s6;
109  cartesian.at(2) = -deviatoric.at(1) / s6 + deviatoric.at(2) / s2;
110  cartesian.at(3) = -deviatoric.at(1) / s6 - deviatoric.at(2) / s2;
111  //
112  cartesian.at(4) = ( deviatoric.at(3) + deviatoric.at(6) ) * 0.5;
113  cartesian.at(5) = ( deviatoric.at(4) + deviatoric.at(7) ) * 0.5;
114  cartesian.at(6) = ( deviatoric.at(5) + deviatoric.at(8) ) * 0.5;
115 }
116 
117 
118 void MixedGradientPressureNeumann :: fromDeviatoricBase2D(FloatMatrix &cartesian, FloatMatrix &deviatoric)
119 {
120  cartesian.resize(3, 3);
121  // E1 = [1/sqrt(2), 1/sqrt(2), 0,0]'; E2 = [0,0,1,0]'; E3 = [0,0,0,1]';
122  // C = E1*(E1'*e11 + E2'*e12 + E3'*e13) + E2*(E1'*e21 + E2'*e22 + E3'*e23) + E3*(E1'*e31 + E2'*e32 + E3'*e33);
123 
124  cartesian.at(1, 1) = deviatoric.at(1, 1) * 0.5;
125  cartesian.at(2, 2) = deviatoric.at(1, 1) * 0.5;
126  cartesian.at(1, 2) = -deviatoric.at(1, 1) * 0.5;
127  cartesian.at(2, 1) = -deviatoric.at(1, 1) * 0.5;
128  //
129  cartesian.at(1, 3) = ( deviatoric.at(1, 2) + deviatoric.at(1, 3) ) / sqrt(8.0);
130  cartesian.at(2, 3) = -( deviatoric.at(1, 2) + deviatoric.at(1, 3) ) / sqrt(8.0);
131  cartesian.at(3, 1) = ( deviatoric.at(2, 1) + deviatoric.at(3, 1) ) / sqrt(8.0);
132  cartesian.at(3, 2) = -( deviatoric.at(2, 1) + deviatoric.at(3, 1) ) / sqrt(8.0);
133  //
134  cartesian.at(3, 3) = ( deviatoric.at(2, 2) + deviatoric.at(2, 3) +
135  deviatoric.at(3, 2) + deviatoric.at(3, 3) ) * 0.25;
136 }
137 
138 
139 void MixedGradientPressureNeumann :: fromDeviatoricBase3D(FloatMatrix &cartesian, FloatMatrix &deviatoric)
140 {
141  cartesian.resize(6, 6);
142  /*
143  * syms
144  * e11 e12 e13 e14 e15 e16 e17 e18
145  * e21 e22 e23 e24 e25 e26 e27 e28
146  * e31 e32 e33 e34 e35 e36 e37 e38
147  * e41 e42 e43 e44 e45 e46 e47 e48
148  * e51 e52 e53 e54 e55 e56 e57 e58
149  * e61 e62 e63 e64 e65 e66 e67 e68
150  * e71 e72 e73 e74 e75 e76 e77 e78
151  * e81 e82 e83 e84 e85 e86 e87 e88;
152  * C = ...
153  * E1*(E1'*e11 + E2'*e12 + E3'*e13 + E4'*e14 + E5'*e15 + E6'*e16 + E7'*e17 + E8'*e18) + ...
154  * E2*(E1'*e21 + E2'*e22 + E3'*e23 + E4'*e24 + E5'*e25 + E6'*e26 + E7'*e27 + E8'*e28) + ...
155  * E3*(E1'*e31 + E2'*e32 + E3'*e33 + E4'*e34 + E5'*e35 + E6'*e36 + E7'*e37 + E8'*e38) + ...
156  * E4*(E1'*e41 + E2'*e42 + E3'*e43 + E4'*e44 + E5'*e45 + E6'*e46 + E7'*e47 + E8'*e48) + ...
157  * E5*(E1'*e51 + E2'*e52 + E3'*e53 + E4'*e54 + E5'*e55 + E6'*e56 + E7'*e57 + E8'*e58) + ...
158  * E6*(E1'*e61 + E2'*e62 + E3'*e63 + E4'*e64 + E5'*e65 + E6'*e66 + E7'*e67 + E8'*e68) + ...
159  * E7*(E1'*e71 + E2'*e72 + E3'*e73 + E4'*e74 + E5'*e75 + E6'*e76 + E7'*e77 + E8'*e78) + ...
160  * E8*(E1'*e81 + E2'*e82 + E3'*e83 + E4'*e84 + E5'*e85 + E6'*e86 + E7'*e87 + E8'*e88);
161  */
162 
163  cartesian.at(1, 1) = deviatoric.at(1, 1) * 2.0 / 3.0;
164  cartesian.at(1, 2) = -deviatoric.at(1, 1) / 3.0 + deviatoric.at(1, 2) * sqrt(3.0);
165  cartesian.at(2, 1) = -deviatoric.at(1, 1) / 3.0 + deviatoric.at(1, 2) * sqrt(3.0);
166  cartesian.at(1, 3) = -deviatoric.at(1, 1) / 3.0 - deviatoric.at(1, 2) * sqrt(3.0);
167  cartesian.at(3, 1) = -deviatoric.at(1, 1) / 3.0 - deviatoric.at(1, 2) * sqrt(3.0);
168 
169  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);
170  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);
171  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);
172  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);
173  // upper off diagonal part
174  cartesian.at(1, 4) = ( deviatoric.at(1, 3) + deviatoric.at(1, 6) ) / sqrt(6.0);
175  cartesian.at(1, 5) = ( deviatoric.at(1, 4) + deviatoric.at(1, 7) ) / sqrt(6.0);
176  cartesian.at(1, 6) = ( deviatoric.at(1, 5) + deviatoric.at(1, 8) ) / sqrt(6.0);
177  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);
178  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);
179  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);
180  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);
181  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);
182  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);
183  // lower off diagonal part
184  cartesian.at(1, 4) = ( deviatoric.at(3, 1) + deviatoric.at(6, 1) ) / sqrt(6.0);
185  cartesian.at(1, 5) = ( deviatoric.at(4, 1) + deviatoric.at(7, 1) ) / sqrt(6.0);
186  cartesian.at(1, 6) = ( deviatoric.at(5, 1) + deviatoric.at(8, 1) ) / sqrt(6.0);
187  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);
188  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);
189  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);
190  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);
191  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);
192  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);
193  //
194  cartesian.at(4, 4) = ( deviatoric.at(3, 3) + deviatoric.at(3, 6) + deviatoric.at(6, 3) + deviatoric.at(6, 6) ) * 0.25;
195  cartesian.at(4, 5) = ( deviatoric.at(3, 4) + deviatoric.at(3, 7) + deviatoric.at(6, 4) + deviatoric.at(6, 7) ) * 0.25;
196  cartesian.at(4, 6) = ( deviatoric.at(3, 5) + deviatoric.at(3, 8) + deviatoric.at(6, 5) + deviatoric.at(6, 8) ) * 0.25;
197  cartesian.at(5, 4) = ( deviatoric.at(4, 3) + deviatoric.at(4, 6) + deviatoric.at(7, 3) + deviatoric.at(7, 6) ) * 0.25;
198  cartesian.at(5, 5) = ( deviatoric.at(4, 4) + deviatoric.at(4, 7) + deviatoric.at(7, 4) + deviatoric.at(7, 7) ) * 0.25;
199  cartesian.at(5, 6) = ( deviatoric.at(4, 5) + deviatoric.at(4, 8) + deviatoric.at(7, 5) + deviatoric.at(7, 8) ) * 0.25;
200  cartesian.at(6, 4) = ( deviatoric.at(5, 3) + deviatoric.at(5, 6) + deviatoric.at(8, 3) + deviatoric.at(8, 6) ) * 0.25;
201  cartesian.at(6, 5) = ( deviatoric.at(5, 4) + deviatoric.at(5, 7) + deviatoric.at(8, 4) + deviatoric.at(8, 7) ) * 0.25;
202  cartesian.at(6, 6) = ( deviatoric.at(5, 5) + deviatoric.at(5, 8) + deviatoric.at(8, 5) + deviatoric.at(8, 8) ) * 0.25;
203 }
204 
205 
206 void MixedGradientPressureNeumann :: setPrescribedDeviatoricGradientFromVoigt(const FloatArray &t)
207 {
208  // Converts the Voigt vector to a deviatoric base dyads (which don't assume symmetry)
209  int nsd = this->domain->giveNumberOfSpatialDimensions();
210  if ( nsd == 3 ) {
211  this->devGradient.resize(8);
212  this->devGradient.at(1) = ( 2.0 * t.at(1) - t.at(2) - t.at(3) ) / sqrt(6.);
213  this->devGradient.at(2) = ( t.at(2) - t.at(3) ) / sqrt(2.);
214  //
215  this->devGradient.at(3) = 0.5 * t.at(4);
216  this->devGradient.at(4) = 0.5 * t.at(5);
217  this->devGradient.at(5) = 0.5 * t.at(6);
218  //
219  this->devGradient.at(6) = 0.5 * t.at(4);
220  this->devGradient.at(7) = 0.5 * t.at(5);
221  this->devGradient.at(8) = 0.5 * t.at(6);
222 
223  this->volGradient = t.at(1) + t.at(2) + t.at(3);
224  } else if ( nsd == 2 ) {
225  this->devGradient.resize(3);
226  this->devGradient.at(1) = ( t.at(1) - t.at(2) ) / sqrt(2.);
227  this->devGradient.at(2) = 0.5 * t.at(3);
228  this->devGradient.at(3) = 0.5 * t.at(3);
229 
230  this->volGradient = t.at(1) + t.at(2);
231  } else {
232  this->devGradient.clear();
233  this->volGradient = t.at(1);
234  }
235 }
236 
237 
238 void MixedGradientPressureNeumann :: giveLocationArrays(std :: vector< IntArray > &rows, std :: vector< IntArray > &cols, CharType type,
239  const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
240 {
241  if ( type != TangentStiffnessMatrix ) {
242  return;
243  }
244 
245  IntArray bNodes;
246  IntArray loc_r, loc_c, sigma_loc_r, sigma_loc_c;
247 
248  // Fetch the columns/rows for the stress contributions;
249  this->sigmaDev->giveLocationArray(dev_id, sigma_loc_r, r_s);
250  this->sigmaDev->giveLocationArray(dev_id, sigma_loc_c, c_s);
251 
252  Set *set = this->giveDomain()->giveSet(this->set);
253  const IntArray &boundaries = set->giveBoundaryList();
254 
255  rows.resize( boundaries.giveSize() );
256  cols.resize( boundaries.giveSize() );
257  int i = 0;
258  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
259  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
260  int boundary = boundaries.at(pos * 2);
261 
262  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
263  e->giveBoundaryLocationArray(loc_r, bNodes, this->dofs, r_s);
264  e->giveBoundaryLocationArray(loc_c, bNodes, this->dofs, c_s);
265  // For most uses, loc_r == loc_c, and sigma_loc_r == sigma_loc_c.
266  rows [ i ] = loc_r;
267  cols [ i ] = sigma_loc_c;
268  i++;
269  // and the symmetric part (usually the transpose of above)
270  rows [ i ] = sigma_loc_r;
271  cols [ i ] = loc_c;
272  i++;
273  }
274 }
275 
276 
277 void MixedGradientPressureNeumann :: integrateVolTangent(FloatArray &answer, Element *e, int boundary)
278 {
279  FloatArray normal, n;
280  FloatMatrix nMatrix;
281 
282  FEInterpolation *interp = e->giveInterpolation(); // Geometry interpolation
283  // Some assumptions here, either velocity or displacement unknowns. Perhaps its better to just as for a certain equation id?
284  // Code will be written using v for velocity, but it could represent displacements.
285  FEInterpolation *interpUnknown = e->giveInterpolation(V_u);
286  if ( interpUnknown ) {
287  interpUnknown = e->giveInterpolation(D_u);
288  }
289 
290  int nsd = e->giveDomain()->giveNumberOfSpatialDimensions();
291  // Order here should be the normal (which takes the first derivative) thus -1
292  int order = interp->giveInterpolationOrder() - 1 + interpUnknown->giveInterpolationOrder();
293  std :: unique_ptr< IntegrationRule > ir( interp->giveBoundaryIntegrationRule(order, boundary) );
294 
295  answer.clear();
296  for ( GaussPoint *gp: *ir ) {
297  const FloatArray &lcoords = gp->giveNaturalCoordinates();
298  FEIElementGeometryWrapper cellgeo(e);
299 
300  // Evaluate the normal;
301  double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
302  // Evaluate the velocity/displacement coefficients
303  interpUnknown->boundaryEvalN(n, boundary, lcoords, cellgeo);
304  nMatrix.beNMatrixOf(n, nsd);
305 
306  answer.plusProduct( nMatrix, normal, detJ * gp->giveWeight() );
307  }
308 }
309 
310 
311 void MixedGradientPressureNeumann :: integrateDevTangent(FloatMatrix &answer, Element *e, int boundary)
312 {
313  FloatArray normal, n;
314  FloatMatrix nMatrix, E_n;
315  FloatMatrix contrib;
316 
317  FEInterpolation *interp = e->giveInterpolation(); // Geometry interpolation
318  // Some assumptions here, either velocity or displacement unknowns. Perhaps its better to just as for a certain equation id?
319  // Code will be written using v for velocity, but it could represent displacements.
320  FEInterpolation *interpUnknown = e->giveInterpolation(V_u);
321  if ( interpUnknown ) {
322  interpUnknown = e->giveInterpolation(D_u);
323  }
324 
325  int nsd = e->giveDomain()->giveNumberOfSpatialDimensions();
326  // Order here should be the normal (which takes the first derivative) thus -1
327  int order = interp->giveInterpolationOrder() - 1 + interpUnknown->giveInterpolationOrder();
328  std :: unique_ptr< IntegrationRule > ir( interp->giveBoundaryIntegrationRule(order, boundary) );
329 
330  answer.clear();
331  for ( auto &gp: *ir ) {
332  const FloatArray &lcoords = gp->giveNaturalCoordinates();
333  FEIElementGeometryWrapper cellgeo(e);
334 
335  // Evaluate the normal;
336  double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
337  // Evaluate the velocity/displacement coefficients
338  interpUnknown->boundaryEvalN(n, boundary, lcoords, cellgeo);
339  nMatrix.beNMatrixOf(n, nsd);
340 
341  // Formulating like this to avoid third order tensors, which is hard to express in linear algebra.
342  // v (x) n : E_i = v . ( E_i . n ) = v . E_n
343  if ( nsd == 3 ) {
344  E_n.resize(8, 3);
345  E_n.at(1, 1) = 2.0 * normal.at(1) / sqrt(6.0);
346  E_n.at(1, 2) = -normal.at(2) / sqrt(6.0);
347  E_n.at(1, 3) = -normal.at(3) / sqrt(6.0);
348 
349  E_n.at(2, 1) = 0.;
350  E_n.at(2, 2) = normal.at(2) / sqrt(2.0);
351  E_n.at(2, 3) = -normal.at(3) / sqrt(2.0);
352 
353  E_n.at(3, 1) = normal.at(2);
354  E_n.at(3, 2) = 0.;
355  E_n.at(3, 3) = 0.;
356 
357  E_n.at(4, 1) = normal.at(3);
358  E_n.at(4, 2) = 0.;
359  E_n.at(4, 3) = 0.;
360 
361  E_n.at(5, 1) = 0.;
362  E_n.at(5, 2) = normal.at(3);
363  E_n.at(5, 3) = 0.;
364 
365  E_n.at(6, 1) = 0.;
366  E_n.at(6, 2) = normal.at(1);
367  E_n.at(6, 3) = 0.;
368 
369  E_n.at(7, 1) = 0.;
370  E_n.at(7, 2) = 0.;
371  E_n.at(7, 3) = normal.at(1);
372 
373  E_n.at(8, 1) = 0.;
374  E_n.at(8, 2) = 0.;
375  E_n.at(8, 3) = normal.at(2);
376  } else if ( nsd == 2 ) {
377  E_n.resize(3, 2);
378  E_n.at(1, 1) = normal.at(1) / sqrt(2.0);
379  E_n.at(1, 2) = -normal.at(2) / sqrt(2.0);
380 
381  E_n.at(2, 1) = normal.at(2);
382  E_n.at(2, 2) = 0.;
383 
384  E_n.at(3, 1) = 0.;
385  E_n.at(3, 2) = normal.at(1);
386  } else {
387  E_n.clear();
388  }
389 
390  contrib.beProductOf(E_n, nMatrix);
391 
392  answer.add(detJ * gp->giveWeight(), contrib);
393  }
394 }
395 
396 
397 void MixedGradientPressureNeumann :: assembleVector(FloatArray &answer, TimeStep *tStep,
398  CharType type, ValueModeType mode, const UnknownNumberingScheme &s, FloatArray *eNorms)
399 {
400  Set *set = this->giveDomain()->giveSet(this->set);
401  const IntArray &boundaries = set->giveBoundaryList();
402 
403  IntArray loc, sigma_loc; // For the velocities and stress respectively
404  IntArray masterDofIDs;
405  IntArray bNodes;
406  this->sigmaDev->giveLocationArray(dev_id, sigma_loc, s);
407 
408  if ( type == ExternalForcesVector ) {
409  // The external forces have two contributions. On the additional equations for sigmaDev, the load is simple the deviatoric gradient.
410  double rve_size = this->domainSize();
411  FloatArray devLoad;
412  devLoad.beScaled(-rve_size, this->devGradient);
413  answer.assemble(devLoad, sigma_loc);
414 
415  // The second contribution is on the momentumbalance equation; - int delta_v . n dA * p
416  FloatArray fe;
417  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
418  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
419  int boundary = boundaries.at(pos * 2);
420 
421  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
422  e->giveBoundaryLocationArray(loc, bNodes, this->dofs, s, & masterDofIDs);
423  this->integrateVolTangent(fe, e, boundary);
424  fe.times(-this->pressure);
425  answer.assemble(fe, loc);
426  if ( eNorms ) {
427  eNorms->assembleSquared(fe, masterDofIDs);
428  }
429  }
430  } else if ( type == InternalForcesVector ) {
431  FloatMatrix Ke;
432  FloatArray fe_v, fe_s;
433  FloatArray s_dev, e_v;
434 
435  // Fetch the current values of the stress;
436  this->sigmaDev->giveUnknownVector(s_dev, dev_id, mode, tStep);
437 
438  // Assemble: int delta_v (x) n dA : E_i s_i
439  // int v (x) n dA : E_i delta_s_i
440  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
441  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
442  int boundary = boundaries.at(pos * 2);
443 
444  // Fetch the element information;
445  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
446  e->giveBoundaryLocationArray(loc, bNodes, this->dofs, s, & masterDofIDs);
447  e->computeBoundaryVectorOf(bNodes, this->dofs, mode, tStep, e_v);
448  this->integrateDevTangent(Ke, e, boundary);
449 
450  // We just use the tangent, less duplicated code (the addition of sigmaDev is linear).
451  fe_v.beProductOf(Ke, e_v);
452  fe_s.beTProductOf(Ke, s_dev);
453  // Note: The terms appear negative in the equations:
454  fe_v.negated();
455  fe_s.negated();
456 
457  answer.assemble(fe_s, loc); // Contributions to delta_v equations
458  answer.assemble(fe_v, sigma_loc); // Contribution to delta_s_i equations
459  if ( eNorms ) {
460  eNorms->assembleSquared(fe_s, masterDofIDs);
461  eNorms->assembleSquared(fe_v, dev_id);
462  }
463  }
464  }
465 }
466 
467 
468 void MixedGradientPressureNeumann :: assemble(SparseMtrx &answer, TimeStep *tStep,
469  CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, double scale)
470 {
471  if ( type == TangentStiffnessMatrix || type == SecantStiffnessMatrix || type == ElasticStiffnessMatrix ) {
472  FloatMatrix Ke, KeT;
473  IntArray loc_r, loc_c, sigma_loc_r, sigma_loc_c;
474  IntArray bNodes;
475  Set *set = this->giveDomain()->giveSet(this->set);
476  const IntArray &boundaries = set->giveBoundaryList();
477 
478  // Fetch the columns/rows for the stress contributions;
479  this->sigmaDev->giveLocationArray(dev_id, sigma_loc_r, r_s);
480  this->sigmaDev->giveLocationArray(dev_id, sigma_loc_c, c_s);
481 
482  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
483  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
484  int boundary = boundaries.at(pos * 2);
485 
486  // Fetch the element information;
487  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
488  e->giveBoundaryLocationArray(loc_r, bNodes, this->dofs, r_s);
489  e->giveBoundaryLocationArray(loc_c, bNodes, this->dofs, c_s);
490  this->integrateDevTangent(Ke, e, boundary);
491  Ke.negated();
492  Ke.times(scale);
493  KeT.beTranspositionOf(Ke);
494  answer.assemble(sigma_loc_r, loc_c, Ke); // Contribution to delta_s_i equations
495  answer.assemble(loc_r, sigma_loc_c, KeT); // Contributions to delta_v equations
496  }
497  }
498 }
499 
500 
501 void MixedGradientPressureNeumann :: computeFields(FloatArray &sigmaDev, double &vol, TimeStep *tStep)
502 {
503  int nsd = this->giveDomain()->giveNumberOfSpatialDimensions();
504  FloatArray sigmaDevBase;
505  Set *set = this->giveDomain()->giveSet(this->set);
506  const IntArray &boundaries = set->giveBoundaryList();
507 
508  // Fetch the current values of the stress in deviatoric base;
509  sigmaDevBase.resize( this->sigmaDev->giveNumberOfDofs() );
510  for ( int i = 1; i <= this->sigmaDev->giveNumberOfDofs(); i++ ) {
511  sigmaDevBase.at(i) = this->sigmaDev->giveDofWithID(dev_id.at(i))->giveUnknown(VM_Total, tStep);
512  }
513  // Convert it back from deviatoric base:
514  if ( nsd == 3 ) {
515  this->fromDeviatoricBase3D(sigmaDev, sigmaDevBase);
516  } else if ( nsd == 2 ) {
517  this->fromDeviatoricBase2D(sigmaDev, sigmaDevBase);
518  } else {
519  sigmaDev.clear();
520  }
521 
522  // Postprocessing; vol = int v . n dA
523  IntArray bNodes;
524  FloatArray unknowns, fe;
525  vol = 0.;
526  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
527  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
528  int boundary = boundaries.at(pos * 2);
529 
530  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
531  e->computeBoundaryVectorOf(bNodes, this->dofs, VM_Total, tStep, unknowns);
532  this->integrateVolTangent(fe, e, boundary);
533  vol += fe.dotProduct(unknowns);
534  }
535  double rve_size = this->domainSize();
536  vol /= rve_size;
537  vol -= volGradient; // This is needed for consistency; We return the volumetric "residual" if a gradient with volumetric contribution is set.
538 }
539 
540 
541 void MixedGradientPressureNeumann :: computeTangents(FloatMatrix &Ed, FloatArray &Ep, FloatArray &Cd, double &Cp, TimeStep *tStep)
542 {
543  // Fetch some information from the engineering model
544  EngngModel *rve = this->giveDomain()->giveEngngModel();
546  std :: unique_ptr< SparseLinearSystemNM > solver(
547  classFactory.createSparseLinSolver( ST_Petsc, this->domain, this->domain->giveEngngModel() ) ); // = rve->giveLinearSolver();
548  SparseMtrxType stype = solver->giveRecommendedMatrix(true);
549  EModelDefaultEquationNumbering fnum;
550  Set *set = this->giveDomain()->giveSet(this->set);
551  IntArray bNodes;
552  const IntArray &boundaries = set->giveBoundaryList();
553  double rve_size = this->domainSize();
554 
555  // Set up and assemble tangent FE-matrix which will make up the sensitivity analysis for the macroscopic material tangent.
556  std :: unique_ptr< SparseMtrx > Kff( classFactory.createSparseMtrx(stype) );
557  if ( !Kff ) {
558  OOFEM_ERROR("Couldn't create sparse matrix of type %d\n", stype);
559  }
560  Kff->buildInternalStructure(rve, this->domain->giveNumber(), fnum);
561  rve->assemble(*Kff, tStep, TangentAssembler(TangentStiffness), fnum, this->domain);
562 
563  // Setup up indices and locations
564  int neq = Kff->giveNumberOfRows();
565 
566  // Indices and such of internal dofs
567  int ndev = this->sigmaDev->giveNumberOfDofs();
568 
569  // Matrices and arrays for sensitivities
570  FloatMatrix ddev_pert(neq, ndev); // In fact, npeq should most likely equal ndev
571  FloatArray p_pert(neq); // RHS for d_dev [d_dev11, d_dev22, d_dev12] in 2D
572 
573  FloatMatrix s_d(neq, ndev); // Sensitivity fields for d_dev
574  FloatArray s_p(neq); // Sensitivity fields for p
575 
576  // Unit pertubations for d_dev
577  ddev_pert.zero();
578  for ( int i = 1; i <= ndev; ++i ) {
579  int eqn = this->sigmaDev->giveDofWithID(dev_id.at(i))->giveEquationNumber(fnum);
580  ddev_pert.at(eqn, i) = -1.0 * rve_size;
581  }
582 
583  // Unit pertubation for d_p
584  p_pert.zero();
585  FloatArray fe;
586  IntArray loc;
587  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
588  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
589  int boundary = boundaries.at(pos * 2);
590 
591  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
592  e->giveBoundaryLocationArray(loc, bNodes, this->dofs, fnum);
593  this->integrateVolTangent(fe, e, boundary);
594  fe.times(-1.0); // here d_p = 1.0
595  p_pert.assemble(fe, loc);
596  }
597 
598  // Solve all sensitivities
599  solver->solve(*Kff, ddev_pert, s_d);
600  solver->solve(*Kff, p_pert, s_p);
601 
602  // Extract the stress response from the solutions
603  FloatArray sigma_p(ndev);
604  FloatMatrix sigma_d(ndev, ndev);
605  for ( int i = 1; i <= ndev; ++i ) {
606  int eqn = this->sigmaDev->giveDofWithID(dev_id.at(i))->giveEquationNumber(fnum);
607  sigma_p.at(i) = s_p.at(eqn);
608  for ( int j = 1; j <= ndev; ++j ) {
609  sigma_d.at(i, j) = s_d.at(eqn, j);
610  }
611  }
612 
613  // Post-process the volumetric rate of deformations in the sensitivity fields;
614  FloatArray e_d(ndev);
615  e_d.zero();
616  double e_p = 0.0;
617  for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
618  Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
619  int boundary = boundaries.at(pos * 2);
620 
621  this->integrateVolTangent(fe, e, boundary);
622 
623  e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);
624  e->giveBoundaryLocationArray(loc, bNodes, this->dofs, fnum);
625 
626  // Using "loc" to pick out the relevant contributions. This won't work at all if there are local coordinate systems in these nodes
627  // or slave nodes etc. The goal is to compute the velocity from the sensitivity field, but we need to avoid going through the actual
628  // engineering model. If this ever becomes an issue it needs to perform the same steps as Element::giveUnknownVector does.
629  for ( int i = 1; i <= fe.giveSize(); ++i ) {
630  if ( loc.at(i) > 0 ) {
631  e_p += fe.at(i) * s_p.at( loc.at(i) );
632  for ( int j = 1; j <= ndev; ++j ) {
633  e_d.at(j) += fe.at(i) * s_d.at(loc.at(i), j);
634  }
635  }
636  }
637  }
638  e_p /= rve_size;
639  e_d.times(1. / rve_size);
640 
641  // Now we need to express the tangents in the normal cartesian coordinate system (as opposed to the deviatoric base we use during computations
642  Cp = - e_p; // Scalar components are of course the same
643  double nsd = this->giveDomain()->giveNumberOfSpatialDimensions();
644  if ( nsd == 3 ) {
645  this->fromDeviatoricBase3D(Cd, e_d);
646  this->fromDeviatoricBase3D(Ep, sigma_p);
647  this->fromDeviatoricBase3D(Ed, sigma_d);
648  } else if ( nsd == 2 ) {
649  this->fromDeviatoricBase2D(Cd, e_d);
650  this->fromDeviatoricBase2D(Ep, sigma_p);
651  this->fromDeviatoricBase2D(Ed, sigma_d);
652  } else { // For 1D case, there simply are no deviatoric components!
653  Cd.clear();
654  Ep.clear();
655  Ed.clear();
656  }
657 }
658 
659 
660 IRResultType MixedGradientPressureNeumann :: initializeFrom(InputRecord *ir)
661 {
662  return MixedGradientPressureBC :: initializeFrom(ir);
663 }
664 
665 
666 void MixedGradientPressureNeumann :: giveInputRecord(DynamicInputRecord &input)
667 {
668  MixedGradientPressureBC :: giveInputRecord(input);
669  input.setField(this->pressure, _IFT_MixedGradientPressure_pressure);
670  OOFEM_ERROR("Not supported yet");
671  //FloatArray devGradientVoigt;
672  //input.setField(devGradientVoigt, _IFT_MixedGradientPressure_devGradient);
673 }
674 
675 
676 
677 void MixedGradientPressureNeumann :: scale(double s)
678 {
679  devGradient.times(s);
680  pressure *= s;
681 }
682 } // end namespace oofem
oofem::MixedGradientPressureNeumann::MixedGradientPressureNeumann
MixedGradientPressureNeumann(int n, Domain *d)
Creates boundary condition with given number, belonging to given domain.
Definition: mixedgradientpressureneumann.C:62
oofem::DynamicInputRecord::setField
void setField(int item, InputFieldType id)
Definition: dynamicinputrecord.C:299
oofem::EModelDefaultEquationNumbering
The representation of EngngModel default unknown numbering.
Definition: unknownnumberingscheme.h:80
oofem::REGISTER_BoundaryCondition
REGISTER_BoundaryCondition(BoundaryCondition)
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:94
oofem::MixedGradientPressureNeumann::scale
virtual void scale(double s)
Scales the receiver according to given value.
Definition: mixedgradientpressureneumann.C:677
oofem::Domain
Class and object Domain.
Definition: domain.h:115
oofem::SparseMtrx::assemble
virtual int assemble(const IntArray &loc, const FloatMatrix &mat)=0
Assembles sparse matrix from contribution of local elements.
oofem::FEMComponent::domain
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition: femcmpnn.h:82
oofem::SparseMtrx
Base class for all matrices stored in sparse format.
Definition: sparsemtrx.h:60
oofem::ClassFactory::createSparseMtrx
SparseMtrx * createSparseMtrx(SparseMtrxType type)
Creates new instance of sparse matrix corresponding to given keyword.
Definition: classfactory.C:105
oofem::MixedGradientPressureNeumann::giveInputRecord
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: mixedgradientpressureneumann.C:666
oofem::FEInterpolation::boundaryEvalNormal
virtual double boundaryEvalNormal(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates the normal on the requested boundary.
oofem::FloatArray::at
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
oofem::ValueModeType
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
oofem::MixedGradientPressureBC::domainSize
double domainSize()
Computes the size (including pores) by surface integral over the domain.
Definition: mixedgradientpressurebc.C:43
oofem::FloatArray::clear
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
oofem::FEInterpolation::giveInterpolationOrder
int giveInterpolationOrder()
Returns the interpolation order.
Definition: feinterpol.h:159
oofem::Domain::giveEngngModel
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
oofem::TangentAssembler
Implementation for assembling tangent matrices in standard monolithic FE-problems.
Definition: assemblercallback.h:185
oofem::Element
Abstract base class for all finite elements.
Definition: element.h:145
oofem::DofManager
Base class for dof managers.
Definition: dofmanager.h:113
oofem::FloatMatrix::negated
void negated()
Changes sign of receiver values.
Definition: floatmatrix.C:1605
oofem::Domain::giveNumber
int giveNumber()
Returns domain number.
Definition: domain.h:266
oofem::ActiveBoundaryCondition::giveUnknown
virtual double giveUnknown(PrimaryField &field, ValueModeType mode, TimeStep *tStep, ActiveDof *dof)
Computes the value of the dof.
Definition: activebc.h:200
oofem::Domain::giveNumberOfSpatialDimensions
int giveNumberOfSpatialDimensions()
Returns number of spatial dimensions.
Definition: domain.C:1067
oofem::IntArray
Class implementing an array of integers.
Definition: intarray.h:61
oofem::IntArray::at
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
oofem::MixedGradientPressureNeumann::giveInternalDofManager
virtual DofManager * giveInternalDofManager(int i)
Returns the volumetric DOF manager for i == 1, and the deviatoric manager for i == 2...
Definition: mixedgradientpressureneumann.C:88
oofem::MixedGradientPressureNeumann::assembleVector
virtual void assembleVector(FloatArray &answer, TimeStep *tStep, CharType type, ValueModeType mode, const UnknownNumberingScheme &s, FloatArray *eNorm=NULL)
Assembles B.C.
Definition: mixedgradientpressureneumann.C:397
oofem::Element::giveInterpolation
virtual FEInterpolation * giveInterpolation() const
Definition: element.h:629
oofem::MixedGradientPressureNeumann::computeTangents
virtual void computeTangents(FloatMatrix &Ed, FloatArray &Ep, FloatArray &Cd, double &Cp, TimeStep *tStep)
Computes the macroscopic tangents through sensitivity analysis.
Definition: mixedgradientpressureneumann.C:541
oofem::FEInterpolation
Class representing a general abstraction for finite element interpolation class.
Definition: feinterpol.h:132
oofem::EngngModel::assemble
virtual void assemble(SparseMtrx &answer, TimeStep *tStep, const MatrixAssembler &ma, const UnknownNumberingScheme &s, Domain *domain)
Assembles characteristic matrix of required type into given sparse matrix.
Definition: engngm.C:803
oofem::MasterDof
Class representing "master" degree of freedom.
Definition: masterdof.h:92
oofem::FloatArray::beScaled
void beScaled(double s, const FloatArray &b)
Sets receiver to be .
Definition: floatarray.C:146
oofem::MixedGradientPressureNeumann::computeFields
virtual void computeFields(FloatArray &sigmaDev, double &vol, TimeStep *tStep)
Computes the homogenized fields through sensitivity analysis.
Definition: mixedgradientpressureneumann.C:501
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:277
oofem::Domain::giveElement
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
oofem::FloatArray::dotProduct
double dotProduct(const FloatArray &x) const
Computes the dot product (or inner product) of receiver and argument.
Definition: floatarray.C:463
oofem::MixedGradientPressureNeumann::volGradient
double volGradient
The volumetric part of what was sent in (needed to return the difference).
Definition: mixedgradientpressureneumann.h:100
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition: error.h:61
oofem::DofIDItem
DofIDItem
Type representing particular dof type.
Definition: dofiditem.h:86
oofem::SparseMtrxType
SparseMtrxType
Enumerative type used to identify the sparse matrix type.
Definition: sparsemtrxtype.h:42
_IFT_MixedGradientPressure_pressure
#define _IFT_MixedGradientPressure_pressure
Definition: mixedgradientpressurebc.h:44
oofem::FloatMatrix::times
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
oofem::Set
Set of elements, boundaries, edges and/or nodes.
Definition: set.h:66
oofem::UnknownNumberingScheme
Abstract base class allowing to control the way, how equations are assigned to individual DOFs...
Definition: unknownnumberingscheme.h:47
oofem::MixedGradientPressureNeumann::devGradient
FloatArray devGradient
Prescribed gradient in Voigt form.
Definition: mixedgradientpressureneumann.h:95
oofem::FEIElementGeometryWrapper
Wrapper around element definition to provide FEICellGeometry interface.
Definition: feinterpol.h:95
oofem::Domain::giveSet
Set * giveSet(int n)
Service for accessing particular domain set.
Definition: domain.C:363
oofem::FloatArray::beProductOf
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
oofem::GeneralBoundaryCondition::dofs
IntArray dofs
Dofs that b.c. is applied to (relevant for Dirichlet type b.c.s).
Definition: generalboundarycondition.h:89
oofem::FloatArray::beTProductOf
void beTProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix^T and anArray.
Definition: floatarray.C:708
oofem::FloatMatrix::at
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
oofem::IntArray::resize
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
oofem::MixedGradientPressureBC::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: mixedgradientpressurebc.C:61
oofem::MixedGradientPressureNeumann::sigmaDev
std::unique_ptr< Node > sigmaDev
DOF-manager containing the unknown deviatoric stress.
Definition: mixedgradientpressureneumann.h:106
oofem::MixedGradientPressureNeumann::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: mixedgradientpressureneumann.C:660
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:103
oofem::FloatMatrix::beNMatrixOf
void beNMatrixOf(const FloatArray &n, int nsd)
Assigns the receiver to be a repeated diagonal matrix.
Definition: floatmatrix.C:505
oofem::FEInterpolation::boundaryGiveNodes
virtual void boundaryGiveNodes(IntArray &answer, int boundary)=0
Gives the boundary nodes for requested boundary number.
oofem::Domain::giveNextFreeDofID
int giveNextFreeDofID(int increment=1)
Gives the next free dof ID.
Definition: domain.C:1411
oofem::FloatArray
Class representing vector of real numbers.
Definition: floatarray.h:82
oofem::ClassFactory::createSparseLinSolver
SparseLinearSystemNM * createSparseLinSolver(LinSystSolverType st, Domain *d, EngngModel *m)
Creates new instance of SparseLinearSystemNM corresponding to given type.
Definition: classfactory.C:120
oofem::FloatMatrix
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
oofem::MixedGradientPressureNeumann::giveLocationArrays
virtual void giveLocationArrays(std::vector< IntArray > &rows, std::vector< IntArray > &cols, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
Gives a list of location arrays that will be assembled.
Definition: mixedgradientpressureneumann.C:238
oofem::IRResultType
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
oofem::FloatArray::assemble
void assemble(const FloatArray &fe, const IntArray &loc)
Assembles the array fe (typically, the load vector of a finite element) into the receiver, using loc as location array.
Definition: floatarray.C:551
oofem::CharType
CharType
Definition: chartype.h:87
oofem::FloatMatrix::resize
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
oofem::InputRecord
Class representing the general Input Record.
Definition: inputrecord.h:101
oofem::MixedGradientPressureNeumann::giveNumberOfInternalDofManagers
virtual int giveNumberOfInternalDofManagers()
Returns the number of internal DOF managers (=2).
Definition: mixedgradientpressureneumann.C:82
oofem::FEInterpolation::giveBoundaryIntegrationRule
virtual IntegrationRule * giveBoundaryIntegrationRule(int order, int boundary)
Sets up a suitable integration rule for integrating over the requested boundary.
Definition: feinterpol.C:63
oofem::FloatMatrix::add
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
oofem::MixedGradientPressureNeumann::assemble
virtual void assemble(SparseMtrx &answer, TimeStep *tStep, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, double scale=1.0)
Assembles B.C.
Definition: mixedgradientpressureneumann.C:468
oofem::IntArray::followedBy
void followedBy(const IntArray &b, int allocChunk=0)
Appends array b at the end of receiver.
Definition: intarray.C:145
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
oofem::FloatArray::assembleSquared
void assembleSquared(const FloatArray &fe, const IntArray &loc)
Assembles the array fe with each component squared.
Definition: floatarray.C:572
oofem::Element::giveBoundaryLocationArray
virtual void giveBoundaryLocationArray(IntArray &locationArray, const IntArray &bNodes, const UnknownNumberingScheme &s, IntArray *dofIds=NULL)
Returns the location array for the boundary of the element.
Definition: element.C:446
oofem::MixedGradientPressureNeumann::dev_id
IntArray dev_id
Dof IDs for the lagrange multipliers in sigmaDev.
Definition: mixedgradientpressureneumann.h:108
oofem::DynamicInputRecord
Class representing the a dynamic Input Record.
Definition: dynamicinputrecord.h:59
oofem::FloatArray::times
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
oofem::FloatMatrix::beTranspositionOf
void beTranspositionOf(const FloatMatrix &src)
Assigns to the receiver the transposition of parameter.
Definition: floatmatrix.C:323
oofem::classFactory
ClassFactory & classFactory
Definition: classfactory.C:59
new
#define new
Definition: oofeggraphiccontext.h:51
oofem::FloatMatrix::zero
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
oofem::FEMComponent::giveDomain
Domain * giveDomain() const
Definition: femcmpnn.h:100
oofem::GeneralBoundaryCondition::giveInputRecord
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: generalboundarycondition.C:108
oofem::EngngModel
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
oofem::FloatMatrix::beProductOf
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
oofem::IntArray::giveSize
int giveSize() const
Definition: intarray.h:203
oofem::FloatArray::giveSize
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
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:311
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
oofem::FloatMatrix::clear
void clear()
Sets size of receiver to be an empty matrix. It will have zero rows and zero columns size...
Definition: floatmatrix.h:516
oofem::Node
Class implementing node in finite element mesh.
Definition: node.h:87
oofem::FloatArray::negated
void negated()
Switches the sign of every coefficient of receiver.
Definition: floatarray.C:739
oofem::MixedGradientPressureBC
General class for boundary condition that prolongates macroscopic fields to incompressible flow...
Definition: mixedgradientpressurebc.h:58
oofem::Element::computeBoundaryVectorOf
void computeBoundaryVectorOf(const IntArray &bNodes, const IntArray &dofIDMask, ValueModeType u, TimeStep *tStep, FloatArray &answer, bool padding=false)
Boundary version of computeVectorOf.
Definition: element.C:139
oofem::FEInterpolation::boundaryEvalN
virtual void boundaryEvalN(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates the basis functions on the requested boundary.
oofem::GaussPoint
Class representing integration point in finite element program.
Definition: gausspoint.h:93
oofem::TimeStep
Class representing solution step.
Definition: timestep.h:80
oofem::MixedGradientPressureNeumann::setPrescribedDeviatoricGradientFromVoigt
virtual void setPrescribedDeviatoricGradientFromVoigt(const FloatArray &ddev)
Sets the prescribed tensor from the matrix from given Voigt notation.
Definition: mixedgradientpressureneumann.C:206
oofem::FloatArray::resize
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
oofem::FloatArray::plusProduct
void plusProduct(const FloatMatrix &b, const FloatArray &s, double dV)
Adds the product .
Definition: floatarray.C:226
This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:30 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011