OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
supgelement.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 "supgelement.h"
36 #include "domain.h"
37 #include "load.h"
38 #include "gausspoint.h"
39 #include "intarray.h"
40 #include "floatarray.h"
41 #include "floatmatrix.h"
42 #include "fluiddynamicmaterial.h"
43 #include "fluidcrosssection.h"
44 #include "dynamicinputrecord.h"
45 #include "engngm.h"
46 #include "node.h"
47 #include "dof.h"
48 
49 #ifdef __OOFEG
50  #include "oofeggraphiccontext.h"
51  #include "connectivitytable.h"
52 #endif
53 
54 namespace oofem {
56  FMElement(n, aDomain), t_supg(0), t_pspg(0), t_lsic(0)
57 { }
58 
59 
61 // Destructor.
62 { }
63 
66 {
67  IRResultType result; // Required by IR_GIVE_FIELD macro
68 
69 
71  if ( !boundarySides.isEmpty() ) {
73  }
74 
75  return FMElement :: initializeFrom(ir);
76 }
77 
78 
79 void
81 {
83  if ( !boundarySides.isEmpty() ) {
86  }
87 }
88 
89 
90 void
92  CharType type, TimeStep *tStep)
93 //
94 // returns characteristics matrix of receiver according to mtrx
95 //
96 {
97 
98  if ( type == TangentStiffnessMatrix ) {
99  // stokes flow only
100  double dscale = this->giveDomain()->giveEngngModel()->giveVariableScale(VST_Density);
101  double uscale = this->giveDomain()->giveEngngModel()->giveVariableScale(VST_Velocity);
102 
103  IntArray vloc, ploc;
104  FloatMatrix h;
105  int size = this->computeNumberOfDofs();
106  this->giveLocalVelocityDofMap(vloc);
107  this->giveLocalPressureDofMap(ploc);
108  answer.resize(size, size);
109  answer.zero();
110 
111  //this->computeAccelerationTerm_MB(h, tStep);
112  //answer.assemble(h, vloc);
113  this->computeDiffusionDerivativeTerm_MB(h, TangentStiffness, tStep);
114  answer.assemble(h, vloc);
115  this->computePressureTerm_MB(h, tStep);
116  answer.assemble(h, vloc, ploc);
117  //this->computeLSICStabilizationTerm_MB(h, tStep);
118  //h.times( alpha * tStep->giveTimeIncrement() * lscale / ( dscale * uscale * uscale ) );
119  //answer.assemble(h, vloc);
120  this->computeBCLhsTerm_MB(h, tStep);
121  if ( h.isNotEmpty() ) {
122  answer.assemble(h, vloc);
123  }
124 
125  this->computeBCLhsPressureTerm_MB(h, tStep);
126  if ( h.isNotEmpty() ) {
127  answer.assemble(h, vloc, ploc);
128  }
129 
130  // conservation eq part
131  this->computeLinearAdvectionTerm_MC(h, tStep);
132  h.times( 1.0 / ( dscale * uscale ) );
133  answer.assemble(h, ploc, vloc);
134  this->computeBCLhsPressureTerm_MC(h, tStep);
135  if ( h.isNotEmpty() ) {
136  answer.assemble(h, ploc, vloc);
137  }
138 
139  this->computeDiffusionDerivativeTerm_MC(h, tStep);
140  answer.assemble(h, ploc, vloc);
141  this->computePressureTerm_MC(h, tStep);
142  answer.assemble(h, ploc);
143  } else {
144  OOFEM_ERROR("giveCharacteristicMatrix: Unknown Type of characteristic mtrx.");
145  }
146 
147 }
148 
149 
150 void
152  TimeStep *tStep)
153 //
154 // returns characteristics vector of receiver according to requested type
155 //
156 {
157  if ( mtrx == ExternalForcesVector ) {
158  answer.clear();
159 #if 0
160  // assembled from assembler from loads directly
161  // stokes flow
162  IntArray vloc, ploc;
163  FloatArray h;
164  int size = this->computeNumberOfDofs();
165  this->giveLocalVelocityDofMap(vloc);
166  this->giveLocalPressureDofMap(ploc);
167  answer.resize(size);
168  answer.zero();
169  this->computeBCRhsTerm_MB(h, tStep);
170  answer.assemble(h, vloc);
171  this->computeBCRhsTerm_MC(h, tStep);
172  answer.assemble(h, ploc);
173 #endif
174  }
175 
176 #if 1
177  else if ( mtrx == InternalForcesVector ) {
178  // stokes flow
179  IntArray vloc, ploc;
180  FloatArray h;
181  int size = this->computeNumberOfDofs();
182  this->giveLocalVelocityDofMap(vloc);
183  this->giveLocalPressureDofMap(ploc);
184  answer.resize(size);
185  answer.zero();
186  //this->computeAdvectionTerm_MB(h, tStep); answer.assemble(h, vloc);
187  this->computeAdvectionTerm_MC(h, tStep);
188  answer.assemble(h, ploc);
189  this->computeDiffusionTerm_MB(h, tStep);
190  answer.assemble(h, vloc);
191  this->computeDiffusionTerm_MC(h, tStep);
192  answer.assemble(h, ploc);
193 
194  FloatMatrix m1;
195  FloatArray v;
196  // add lsic stabilization term
197  //this->giveCharacteristicMatrix(m1, LSICStabilizationTerm_MB, tStep);
198  //m1.times( lscale / ( dscale * uscale * uscale ) );
199  this->computeVectorOfVelocities(VM_Total, tStep, v);
200  //h.beProductOf(m1, v);
201  //answer.assemble(h, vloc);
202  this->computeLinearAdvectionTerm_MC(m1, tStep);
203  //m1.times( 1. / ( dscale * uscale ) );
204  h.beProductOf(m1, v);
205  answer.assemble(h, ploc);
206 
207  // add pressure term
208  this->computePressureTerm_MB(m1, tStep);
209  this->computeVectorOfPressures(VM_Total, tStep, v);
210  h.beProductOf(m1, v);
211  answer.assemble(h, vloc);
212 
213  // pressure term
214  this->computePressureTerm_MC(m1, tStep);
215  this->computeVectorOfPressures(VM_Total, tStep, v);
216  h.beProductOf(m1, v);
217  answer.assemble(h, ploc);
218  }
219 #endif
220  else {
221  OOFEM_ERROR("Unknown Type of characteristic mtrx.");
222  }
223 }
224 
225 
226 
227 
228 
229 void
231 {
232  FloatArray eps;
233 
234  // compute deviatoric strain
235  this->computeDeviatoricStrain(eps, gp, tStep);
236  // call material to compute stress
237  static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->computeDeviatoricStressVector(answer, gp, eps, tStep);
238 }
239 
240 
241 
242 void
244 {
245  bcType boundarytype;
246  int nLoads = 0;
247  //bcType loadtype;
248  FloatMatrix helpMatrix;
249  // loop over boundary load array
250 
251  answer.clear();
252 
253  nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
254  if ( nLoads ) {
255  for ( int i = 1; i <= nLoads; i++ ) {
256  int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
257  int side = boundaryLoadArray.at(i * 2);
258  Load *load = domain->giveLoad(n);
259  boundarytype = load->giveType();
260  if ( boundarytype == SlipWithFriction ) {
261  this->computeSlipWithFrictionBCTerm_MB(helpMatrix, load, side, tStep);
262  answer.add(helpMatrix);
263  } else if ( boundarytype == PenetrationWithResistance ) {
264  this->computePenetrationWithResistanceBCTerm_MB(helpMatrix, load, side, tStep);
265  answer.add(helpMatrix);
266  } else {
267  // OOFEM_ERROR("unsupported load type class");
268  }
269  }
270  }
271 
272  nLoads = this->giveBodyLoadArray()->giveSize();
273 
274  if ( nLoads ) {
275  bcGeomType ltype;
276  for ( int i = 1; i <= nLoads; i++ ) {
277  Load *load = domain->giveLoad( bodyLoadArray.at(i) );
278  ltype = load->giveBCGeoType();
279  if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ReinforceBVT ) ) {
280  this->computeHomogenizedReinforceTerm_MB(helpMatrix, load, tStep);
281  answer.add(helpMatrix);
282  }
283  }
284  }
285 }
286 
287 void
289 {
290  bcType boundarytype;
291  int nLoads = 0;
292  //bcType loadtype;
293  FloatMatrix helpMatrix;
294  // loop over boundary load array
295  answer.clear();
296 
297  nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
298 
299  if ( nLoads ) {
300  for ( int i = 1; i <= nLoads; i++ ) {
301  int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
302  int side = boundaryLoadArray.at(i * 2);
303  Load *load = domain->giveLoad(n);
304  boundarytype = load->giveType();
305  if ( boundarytype == OutFlowBC ) {
306  this->computeOutFlowBCTerm_MB(helpMatrix, side, tStep);
307  answer.add(helpMatrix);
308  } else {
309  //_warning("computeForceLoadVector : unsupported load type class");
310  }
311  }
312  }
313 }
314 
315 void
317 {
318  int nLoads = 0;
319  //bcType loadtype;
320  FloatMatrix helpMatrix;
321 
322  nLoads = this->giveBodyLoadArray()->giveSize();
323  answer.clear();
324  if ( nLoads ) {
325  bcGeomType ltype;
326  for ( int i = 1; i <= nLoads; i++ ) {
327  Load *load = domain->giveLoad( bodyLoadArray.at(i) );
328  ltype = load->giveBCGeoType();
329  if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ReinforceBVT ) ) {
330  this->computeHomogenizedReinforceTerm_MC(helpMatrix, load, tStep);
331  answer.add(helpMatrix);
332  }
333  }
334  }
335 }
336 
337 
338 int
340 //
341 // check internal consistency
342 // mainly tests, whether material and crossSection data
343 // are safe for conversion to "Structural" versions
344 //
345 {
346  int result = 1;
347  return result;
348 }
349 
350 void
352 {
353  FloatArray stress;
354 
355  // force updating strains & stresses
356  for ( auto &iRule: integrationRulesArray ) {
357  for ( GaussPoint *gp: *iRule ) {
358  computeDeviatoricStress(stress, gp, tStep);
359  }
360  }
361 }
362 
363 
364 #ifdef __OOFEG
365 int
367  int node, TimeStep *tStep)
368 {
369  int indx = 1;
370  Node *n = this->giveNode(node);
371 
372  if ( type == IST_Velocity ) {
373  answer.resize( this->giveSpatialDimension() );
374  std::vector< Dof* >::const_iterator dofindx;
375  if ( ( dofindx = n->findDofWithDofId(V_u) ) != n->end() ) {
376  answer.at(indx++) = (*dofindx)->giveUnknown(VM_Total, tStep);
377  } else if ( ( dofindx = n->findDofWithDofId(V_v) ) != n->end() ) {
378  answer.at(indx++) = (*dofindx)->giveUnknown(VM_Total, tStep);
379  } else if ( ( dofindx = n->findDofWithDofId(V_w) ) != n->end() ) {
380  answer.at(indx++) = (*dofindx)->giveUnknown(VM_Total, tStep);
381  }
382 
383  return 1;
384  } else if ( type == IST_Pressure ) {
385  auto dofindx = n->findDofWithDofId(P_f);
386  if ( dofindx != n->end() ) {
387  answer.resize(1);
388  answer.at(1) = (*dofindx)->giveUnknown(VM_Total, tStep);
389  return 1;
390  } else {
391  return 0;
392  }
393  } else {
394  return Element :: giveInternalStateAtNode(answer, type, mode, node, tStep);
395  }
396 }
397 
398 #endif
399 
400 } // end namespace oofem
CrossSection * giveCrossSection()
Definition: element.C:495
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
void setField(int item, InputFieldType id)
IntArray * giveBoundaryLoadArray()
Returns array containing load numbers of boundary loads acting on element.
Definition: element.C:381
virtual void computePressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)=0
Computes pressure terms for momentum balance equations(s).
virtual void giveCharacteristicMatrix(FloatMatrix &answer, CharType type, TimeStep *tStep)
Computes characteristic matrix of receiver of requested type in given time step.
Definition: supgelement.C:91
Class and object Domain.
Definition: domain.h:115
virtual int computeNumberOfDofs()
Computes or simply returns total number of element&#39;s local DOFs.
Definition: element.h:424
Fluid cross-section.
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition: femcmpnn.h:82
virtual void computeLinearAdvectionTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes the linear advection term for mass conservation equation.
bool isEmpty() const
Checks if receiver is empty (i.e., zero sized).
Definition: intarray.h:208
virtual void giveCharacteristicVector(FloatArray &answer, CharType type, ValueModeType mode, TimeStep *tStep)
Computes characteristic vector of receiver of requested type in given time step.
Definition: supgelement.C:151
bcGeomType
Type representing the geometric character of loading.
Definition: bcgeomtype.h:40
virtual void updateInternalState(TimeStep *tStep)
Updates element state after equilibrium in time step has been reached.
Definition: supgelement.C:351
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
#define _IFT_SUPGElement_bcodes
Definition: supgelement.h:44
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
IntArray boundarySides
Array of boundary sides.
Definition: supgelement.h:62
virtual ~SUPGElement()
Definition: supgelement.C:60
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: element.C:706
virtual void computeHomogenizedReinforceTerm_MC(FloatMatrix &answer, Load *load, TimeStep *tStep)
Definition: supgelement.h:190
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: supgelement.C:80
virtual void giveLocalPressureDofMap(IntArray &map)
Definition: supgelement.h:209
virtual void computeBCRhsTerm_MC(FloatArray &answer, TimeStep *tStep)=0
Computes Rhs terms due to boundary conditions.
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
Distributed body load.
Definition: bcgeomtype.h:43
bool isNotEmpty() const
Tests for empty matrix.
Definition: floatmatrix.h:162
virtual void computeVectorOfPressures(ValueModeType mode, TimeStep *tStep, FloatArray &pressures)
Definition: fmelement.C:55
virtual void computeDiffusionTerm_MB(FloatArray &answer, TimeStep *tStep)=0
Computes diffusion terms for momentum balance equations(s).
virtual void computeHomogenizedReinforceTerm_MB(FloatMatrix &answer, Load *load, TimeStep *tStep)
Definition: supgelement.h:186
bcType
Type representing the type of bc.
Definition: bctype.h:40
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: element.C:638
virtual void computeDiffusionTerm_MC(FloatArray &answer, TimeStep *tStep)=0
Computes diffusion terms for mass conservation equation.
virtual void computeVectorOfVelocities(ValueModeType mode, TimeStep *tStep, FloatArray &velocities)
Definition: fmelement.C:48
virtual void computeBCRhsTerm_MB(FloatArray &answer, TimeStep *tStep)=0
Computes Rhs terms due to boundary conditions.
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual double giveVariableScale(VarScaleType varId)
Returns the scale factor for given variable type.
Definition: engngm.h:1087
virtual void computeDeviatoricStress(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
Definition: supgelement.C:230
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
virtual void computePenetrationWithResistanceBCTerm_MB(FloatMatrix &answer, Load *load, int side, TimeStep *tStep)
Computes Lhs contribution due to applied Penetration bc.
Definition: supgelement.h:173
IntArray * giveBodyLoadArray()
Returns array containing load numbers of loads acting on element.
Definition: element.C:372
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
IntArray bodyLoadArray
Array containing indexes of loads (body loads and boundary loads are kept separately), that apply on receiver.
Definition: element.h:160
virtual int giveInternalStateAtNode(FloatArray &answer, InternalStateType type, InternalStateMode mode, int node, TimeStep *tStep)
Returns internal state variable (like stress,strain) at node of element in Reduced form...
Definition: element.C:1661
virtual void computeDeviatoricStrain(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)=0
virtual void computeDiffusionDerivativeTerm_MB(FloatMatrix &answer, MatResponseMode mode, TimeStep *tStep)=0
Computes the derivative of diffusion terms for momentum balance equations(s) with respect to nodal ve...
virtual void computeOutFlowBCTerm_MB(FloatMatrix &answer, int side, TimeStep *tStep)
Computes Lhs contribution due to outflow BC.
Definition: supgelement.h:180
virtual void giveLocalVelocityDofMap(IntArray &map)
Definition: supgelement.h:208
Class representing vector of real numbers.
Definition: floatarray.h:82
This abstract class represent a general base element class for fluid dynamic problems.
Definition: fmelement.h:54
virtual void computeBCLhsPressureTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - pressure.
Definition: supgelement.C:316
virtual int giveSpatialDimension()
Returns the element spatial dimension (1, 2, or 3).
Definition: element.C:1347
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
virtual void computeDiffusionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes diffusion derivative terms for mass conservation equation.
virtual void computeBCLhsPressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - pressure.
Definition: supgelement.C:288
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
CharType
Definition: chartype.h:87
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
Class representing the general Input Record.
Definition: inputrecord.h:101
std::vector< Dof * >::const_iterator findDofWithDofId(DofIDItem dofID) const
Finds index of DOF with required physical meaning of receiver.
Definition: dofmanager.C:266
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual bcGeomType giveBCGeoType() const
Returns geometry character of boundary condition.
IntArray boundaryCodes
Boundary sides codes.
Definition: supgelement.h:64
Class representing the a dynamic Input Record.
std::vector< std::unique_ptr< IntegrationRule > > integrationRulesArray
List of integration rules of receiver (each integration rule contains associated integration points a...
Definition: element.h:170
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
Domain * giveDomain() const
Definition: femcmpnn.h:100
#define _IFT_SUPGElement_bsides
Definition: supgelement.h:43
Load is base abstract class for all loads.
Definition: load.h:61
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: supgelement.C:65
SUPGElement(int n, Domain *aDomain)
Definition: supgelement.C:55
int giveSize() const
Definition: intarray.h:203
the oofem namespace is to define a context or scope in which all oofem names are defined.
std::vector< Dof * >::iterator end()
Definition: dofmanager.h:158
void assemble(const FloatMatrix &src, const IntArray &loc)
Assembles the contribution using localization array into receiver.
Definition: floatmatrix.C:251
void clear()
Sets size of receiver to be an empty matrix. It will have zero rows and zero columns size...
Definition: floatmatrix.h:516
virtual bcValType giveBCValType() const
Returns receiver load type.
Class implementing node in finite element mesh.
Definition: node.h:87
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
virtual int checkConsistency()
Performs consistency check.
Definition: supgelement.C:339
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
Load * giveLoad(int n)
Service for accessing particular domain load.
Definition: domain.C:222
Class representing integration point in finite element program.
Definition: gausspoint.h:93
IntArray boundaryLoadArray
Definition: element.h:160
int giveInternalStateAtNode(FloatArray &answer, InternalStateType type, InternalStateMode mode, int node, TimeStep *tStep)
Returns internal state variable (like stress,strain) at node of element in Reduced form...
Definition: supgelement.C:366
virtual void computeAdvectionTerm_MC(FloatArray &answer, TimeStep *tStep)=0
Computes advection terms for mass conservation equation.
Class representing solution step.
Definition: timestep.h:80
virtual void computePressureTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes pressure terms for mass conservation equation.
InternalStateMode
Determines the mode of internal variable.
virtual void computeBCLhsTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - velocity.
Definition: supgelement.C:243
virtual void computeSlipWithFrictionBCTerm_MB(FloatMatrix &answer, Load *load, int side, TimeStep *tStep)
Computes Lhs term due to applied slip with friction bc.
Definition: supgelement.h:166
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

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:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011