OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
structuralinterfacematerial.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 
37 #include "dynamicinputrecord.h"
38 #include "gausspoint.h"
39 
40 namespace oofem {
42 {
43  this->useNumericalTangent = false;
44 }
45 
46 
47 int
49 {
51  if ( type == IST_InterfaceJump ) {
52  answer = status->giveJump();
53  answer.resizeWithValues(3); // In case some model is not storing all components.
54  return 1;
55  } else if ( type == IST_InterfaceTraction ) {
56  answer = status->giveTraction();
57  answer.resizeWithValues(3);
58  return 1;
59  } else if ( type == IST_InterfaceFirstPKTraction ) {
60  answer = status->giveFirstPKTraction();
61  answer = status->giveTempFirstPKTraction();
62  answer.resizeWithValues(3);
63  return 1;
64  } else if ( type == IST_DeformationGradientTensor ) {
65  answer.beVectorForm( status->giveF() );
66  return 1;
67  } else if ( type == IST_InterfaceNormal ) {
68  answer = status->giveNormal();
69  return 1;
70  } else {
71  return Material :: giveIPValue(answer, gp, type, tStep);
72  }
73 }
74 
75 
78 {
79  IRResultType result; // Required by IR_GIVE_FIELD macro
80 
82 
83  return IRRT_OK;
84 }
85 
86 
87 void
89 {
91 }
92 
93 
94 #if 0
95 void
96 StructuralInterfaceMaterial :: giveStiffnessMatrix_dTdj_Num(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
97 {
98  // Default implementation for computation of the numerical tangent
99  // Computes the material stiffness using a central difference method
100 
102  if ( status ) {
103  FloatMatrix F;
104  F = status->giveTempF();
105  int dim = F.giveNumberOfRows();
106  answer.resize(dim, dim);
107  answer.zero();
108  const double eps = 1.0e-9;
109  FloatArray T, TPlus, TMinus;
110 
111  FloatArray jump, jumpPlus, jumpMinus, Kcolumn;
112  jump = status->giveTempJump();
113  for ( int i = 1; i <= dim; i++ ) {
114  jumpPlus = jumpMinus = jump;
115  jumpPlus.at(i) += eps;
116  jumpMinus.at(i) -= eps;
117  this->giveFirstPKTraction_3d(TPlus, gp, jumpPlus, F, tStep);
118  this->giveFirstPKTraction_3d(TMinus, gp, jumpMinus, F, tStep);
119 
120  Kcolumn = ( TPlus - TMinus );
121  answer.setColumn(Kcolumn, i);
122  }
123  answer.times( 1.0 / ( 2 * eps ) );
124  this->giveFirstPKTraction_3d(T, gp, jump, F, tStep); // reset temp values by recomputing the stress
125  }
126 }
127 #endif
128 
129 
130 void
132 {
133  FloatMatrix answer3D;
134  give3dStiffnessMatrix_Eng(answer3D, mode, gp, tStep);
135  answer.resize(1, 1);
136  answer.at(1, 1) = answer3D.at(1, 1);
137 }
138 
139 
140 void
142 {
143  FloatMatrix answer3D;
144  give3dStiffnessMatrix_Eng(answer3D, mode, gp, tStep);
145  IntArray mask = {1, 2};
146  answer.beSubMatrixOf(answer3D, mask, mask);
147 
148 #if 0
149  answer3D.printYourself("analytical tangent");
150 
151  FloatMatrix answerNum;
152  give2dStiffnessMatrix_Eng_Num(answerNum, mode, gp, tStep);
153  answerNum.printYourself("numerical tangent");
154 
155  FloatMatrix comp;
156  comp = answer3D;
157  comp.subtract(answerNum);
158  comp.printYourself("difference in numerical tangent to mat method");
159 #endif
160 }
161 
162 void
164 {
165  // Default implementation. Will use large deformation dT/dj as stiffness
166  give3dStiffnessMatrix_dTdj(answer, mode, gp, tStep);
167 }
168 
169 
170 void
172 {
173  this->give3dStiffnessMatrix_dTdj_Num(answer, gp, tStep);
174  answer.resize(1, 1);
175  answer.at(1, 1) = answer.at(1, 1);
176 }
177 
178 
179 void
181 {
182  this->give3dStiffnessMatrix_dTdj_Num(answer, gp, tStep);
183  answer.resize(2, 2);
184 // answer.at(1, 1) = answer.at(1, 1);
185 // answer.at(1, 2) = answer.at(1, 3);
186 // answer.at(2, 1) = answer.at(3, 1);
187 // answer.at(2, 2) = answer.at(3, 3);
188 }
189 
190 
191 void
193 {
194  this->give3dStiffnessMatrix_dTdj_Num(answer, gp, tStep);
195  OOFEM_WARNING("Using numerical tangent");
196 }
197 
198 
199 void
201  const FloatMatrix &reducedF, TimeStep *tStep)
202 {
203  FloatArray traction3D;
204  FloatMatrix F(3, 3);
205  F.at(1, 1) = reducedF.at(1, 1);
206  F.at(2, 2) = 1.;
207  F.at(3, 3) = 1.;
208 
209  this->giveFirstPKTraction_3d(traction3D, gp, { jump.at(1), 0., 0. }, F, tStep);
210  answer = FloatArray{ traction3D.at(1) };
211 }
212 
213 void
215  const FloatMatrix &reducedF, TimeStep *tStep)
216 {
217  FloatArray traction3D;
218  FloatMatrix F(3, 3);
219 // F.at(1, 1) = reducedF.at(1, 1);
220 // F.at(1, 3) = reducedF.at(1, 2);
221 // F.at(2, 2) = 1.;
222 // F.at(3, 1) = reducedF.at(2, 1);
223 // F.at(3, 3) = reducedF.at(2, 2);
224  F.at(1, 1) = reducedF.at(1, 1);
225  F.at(1, 2) = reducedF.at(1, 2);
226  F.at(2, 1) = reducedF.at(2, 1);
227  F.at(2, 2) = reducedF.at(2, 2);
228  F.at(3, 3) = 1.;
229 
230  //this->giveFirstPKTraction_3d(traction3D, gp, { jump.at(1), 0., jump.at(2) }, F, tStep);
231  this->giveFirstPKTraction_3d(traction3D, gp, { jump.at(1), jump.at(2), 0. }, F, tStep);
232  answer = { traction3D.at(1), traction3D.at(2) };
233 }
234 
235 void
237 {
238  FloatArray traction3D;
239  this->giveEngTraction_3d(traction3D, gp, { jump.at(1), 0., 0.}, tStep);
240  answer = FloatArray{ traction3D.at(1) };
241 }
242 
243 void
245 {
246  FloatArray traction3D;
247  //this->giveEngTraction_3d(traction3D, gp, { jump.at(1), 0.0, jump.at(2) }, tStep);
248  this->giveEngTraction_3d(traction3D, gp, { jump.at(1), jump.at(2), 0. }, tStep);
249  answer = { traction3D.at(1), traction3D.at(2) };
250 }
251 
252 void
254 {
255  // Default implementation calls first PK version with F = I
256  FloatMatrix F(3, 3);
257  F.beUnitMatrix();
258  giveFirstPKTraction_3d(answer, gp, jump, F, tStep);
259 }
260 
261 
262 
263 // Numerical tangents
264 void
266 {
267  // Default implementation for computation of the numerical tangent
268  // Computes the material stiffness using a central difference method
269 
270  StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus( gp ) );
271  double eps = 1.0e-9;
272  const FloatMatrix &F = status->giveTempF();
273  FloatArray T, TPlus, TMinus;
274  FloatArray jumpPlus, jumpMinus, Kcolumn;
275  FloatArray jump = FloatArray{status->giveTempJump().at(3)};
276  int dim = jump.giveSize();
277  answer.resize( dim, dim );
278 
279  for ( int i = 1; i <= dim; i++ ) {
280  jumpPlus = jumpMinus = jump;
281  jumpPlus.at( i ) += eps;
282  jumpMinus.at( i ) -= eps;
283  this->giveFirstPKTraction_1d( TPlus, gp, jumpPlus, F, tStep );
284  this->giveFirstPKTraction_1d( TMinus, gp, jumpMinus, F, tStep );
285 
286  Kcolumn.beDifferenceOf(TPlus, TMinus);
287  answer.setColumn( Kcolumn, i );
288  }
289  answer.times( 1.0 / ( 2 * eps ) );
290  this->giveFirstPKTraction_1d(T, gp, jump, F, tStep); // reset temp values by recomputing the stress
291 }
292 
293 void
295 {
296  // Default implementation for computation of the numerical tangent
297  // Computes the material stiffness using a central difference method
298 
299  StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus( gp ) );
300  double eps = 1.0e-9;
301  const FloatMatrix &F = status->giveTempF();
302  FloatArray T, TPlus, TMinus;
303  FloatArray jumpPlus, jumpMinus, Kcolumn;
304  FloatArray jump = {status->giveTempJump().at(1), status->giveTempJump().at(2)};
305  int dim = jump.giveSize();
306  answer.resize( dim, dim );
307 
308  for ( int i = 1; i <= dim; i++ ) {
309  jumpPlus = jumpMinus = jump;
310  jumpPlus.at( i ) += eps;
311  jumpMinus.at( i ) -= eps;
312  this->giveFirstPKTraction_2d(TPlus, gp, jumpPlus, F, tStep);
313  this->giveFirstPKTraction_2d(TMinus, gp, jumpMinus, F, tStep);
314 
315  Kcolumn.beDifferenceOf(TPlus, TMinus);
316  answer.setColumn(Kcolumn, i);
317  }
318  answer.times( 1.0 / ( 2 * eps ) );
319  this->giveFirstPKTraction_2d(T, gp, jump, F, tStep); // reset temp values by recomputing the stress
320 }
321 
322 void
324 {
325  // Default implementation for computation of the numerical tangent
326  // Computes the material stiffness using a central difference method
327 
328  StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus( gp ) );
329  double eps = 1.0e-9;
330  const FloatMatrix &F = status->giveTempF();
331  FloatArray T, TPlus, TMinus;
332  FloatArray jumpPlus, jumpMinus, Kcolumn;
333  const FloatArray &jump = status->giveTempJump();
334  int dim = jump.giveSize();
335  answer.resize( dim, dim );
336 
337  for(int i = 1; i <= dim; i++) {
338  jumpPlus = jumpMinus = jump;
339  jumpPlus.at( i ) += eps;
340  jumpMinus.at( i ) -= eps;
341  this->giveFirstPKTraction_3d( TPlus, gp, jumpPlus, F, tStep );
342  this->giveFirstPKTraction_3d( TMinus, gp, jumpMinus, F, tStep );
343 
344  Kcolumn.beDifferenceOf(TPlus, TMinus);
345  answer.setColumn(Kcolumn, i);
346  }
347  answer.times( 1.0 / ( 2 * eps ) );
348  this->giveFirstPKTraction_3d(T, gp, jump, F, tStep); // reset temp values by recomputing the stress
349 }
350 
351 
352 void
354 {
355  // Default implementation for computation of the numerical tangent d(sig)/d(jump)
356  // Computes the material stiffness using a central difference method
357 
358  StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus( gp ) );
359  double eps = 1.0e-9;
360  FloatArray t, tPlus, tMinus;
361  FloatArray jumpPlus, jumpMinus, Kcolumn;
362  FloatArray jump = FloatArray{status->giveTempJump().at(1)};
363  int dim = jump.giveSize();
364  answer.resize( dim, dim );
365  answer.zero();
366 
367  for(int i = 1; i <= dim; i++) {
368  jumpPlus = jumpMinus = jump;
369  jumpPlus.at( i ) += eps;
370  jumpMinus.at( i ) -= eps;
371  this->giveEngTraction_1d(tPlus, gp, jumpPlus, tStep);
372  this->giveEngTraction_1d(tMinus, gp, jumpMinus, tStep);
373 
374  Kcolumn.beDifferenceOf(tPlus, tMinus);
375  answer.setColumn(Kcolumn, i);
376  }
377  answer.times( 1.0 / ( 2 * eps ) );
378  this->giveEngTraction_1d( t, gp, jump, tStep ); // reset temp values by recomputing the stress
379 }
380 
381 void
383 {
384  // Default implementation for computation of the numerical tangent d(sig)/d(jump)
385  // Computes the material stiffness using a central difference method
386 
387  StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus( gp ) );
388  double eps = 1.0e-12;
389  FloatArray t, tPlus, tMinus;
390  FloatArray tempJump, jumpPlus, jumpMinus, Kcolumn;
391  FloatArray jump = {status->giveTempJump().at(1), status->giveTempJump().at(2)};
392  int dim = jump.giveSize();
393  answer.resize(dim, dim);
394  answer.zero();
395 
396  for(int i = 1; i <= dim; i++) {
397  jumpPlus = jumpMinus = jump;
398  jumpPlus.at( i ) += eps;
399  jumpMinus.at( i ) -= eps;
400  this->giveEngTraction_2d(tPlus, gp, jumpPlus, tStep);
401  this->giveEngTraction_2d(tMinus, gp, jumpMinus, tStep);
402 
403  Kcolumn.beDifferenceOf(tPlus, tMinus);
404  answer.setColumn(Kcolumn, i);
405  }
406  answer.times( 1.0 / ( 2 * eps ) );
407  this->giveEngTraction_2d( t, gp, jump, tStep ); // reset temp values by recomputing the stress
408 }
409 
410 void
412 {
413  // Default implementation for computation of the numerical tangent d(sig)/d(jump)
414  // Computes the material stiffness using a central difference method
415 
416  StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus( gp ) );
417  double eps = 1.0e-9;
418  FloatArray t, tPlus, tMinus;
419  FloatArray jumpPlus, jumpMinus, Kcolumn;
420  const FloatArray &jump = status->giveTempJump();
421  int dim = jump.giveSize();
422  answer.resize(dim, dim);
423  answer.zero();
424 
425  for(int i = 1; i <= dim; i++) {
426  jumpPlus = jumpMinus = jump;
427  jumpPlus.at( i ) += eps;
428  jumpMinus.at( i ) -= eps;
429  this->giveEngTraction_3d(tPlus, gp, jumpPlus, tStep);
430  this->giveEngTraction_3d(tMinus, gp, jumpMinus, tStep);
431 
432  Kcolumn.beDifferenceOf(tPlus, tMinus);
433  answer.setColumn(Kcolumn, i);
434  }
435  answer.times( 1.0 / ( 2 * eps ) );
436  this->giveEngTraction_3d( t, gp, jump, tStep ); // reset temp values by recomputing the stress
437 }
438 
439 } // end namespace oofem
virtual void giveEngTraction_2d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
virtual void giveEngTraction_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
void give3dStiffnessMatrix_dTdj_Num(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
Class and object Domain.
Definition: domain.h:115
virtual void give3dStiffnessMatrix_dTdj(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
virtual void giveEngTraction_1d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
void beSubMatrixOf(const FloatMatrix &src, int topRow, int bottomRow, int topCol, int bottomCol)
Assigns to the receiver the sub-matrix of another matrix.
Definition: floatmatrix.C:962
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
void beVectorForm(const FloatMatrix &aMatrix)
Reciever will be a vector with 9 components formed from a 3x3 matrix.
Definition: floatarray.C:992
const FloatArray & giveFirstPKTraction() const
Returns the const pointer to receiver&#39;s first Piola-Kirchhoff traction vector.
void give1dStiffnessMatrix_dTdj_Num(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: material.C:110
Class implementing an array of integers.
Definition: intarray.h:61
MatResponseMode
Describes the character of characteristic material matrix.
virtual void give1dStiffnessMatrix_dTdj(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
Gives the tangent: .
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
void give2dStiffnessMatrix_dTdj_Num(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
virtual void giveFirstPKTraction_2d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, const FloatMatrix &reducedF, TimeStep *tStep)
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
virtual void giveFirstPKTraction_1d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, const FloatMatrix &reducedF, TimeStep *tStep)
Computes the first Piola-Kirchoff traction vector for given total jump/gap and integration point...
virtual void giveFirstPKTraction_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, const FloatMatrix &F, TimeStep *tStep)
const FloatMatrix & giveTempF() const
Returns the const pointer to receiver&#39;s temporary deformation gradient vector.
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
virtual void give2dStiffnessMatrix_Eng(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
This class implements a structural interface material status information.
void resizeWithValues(int s, int allocChunk=0)
Checks size of receiver towards requested bounds.
Definition: floatarray.C:615
Abstract base class for all material models.
Definition: material.h:95
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
const FloatArray & giveJump() const
Returns the const pointer to receiver&#39;s jump.
void give2dStiffnessMatrix_Eng_Num(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
virtual void give1dStiffnessMatrix_Eng(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
void subtract(const FloatMatrix &a)
Subtracts matrix from the receiver.
Definition: floatmatrix.C:1084
Class representing vector of real numbers.
Definition: floatarray.h:82
const FloatArray & giveTraction() const
Returns the const pointer to receiver&#39;s traction vector.
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
StructuralInterfaceMaterial(int n, Domain *d)
Constructor.
void give3dStiffnessMatrix_Eng_Num(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
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
const FloatArray & giveTempFirstPKTraction() const
Returns the const pointer to receiver&#39;s temporary first Piola-Kirchhoff traction vector.
void setColumn(const FloatArray &src, int c)
Sets the values of the matrix in specified column.
Definition: floatmatrix.C:648
Class representing the a dynamic Input Record.
#define _IFT_StructuralInterfaceMaterial_useNumericalTangent
void printYourself() const
Prints matrix to stdout. Useful for debugging.
Definition: floatmatrix.C:1458
virtual void give2dStiffnessMatrix_dTdj(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
void beUnitMatrix()
Sets receiver to unity matrix.
Definition: floatmatrix.C:1332
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
void give1dStiffnessMatrix_Eng_Num(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual void give3dStiffnessMatrix_Eng(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
const FloatArray & giveTempJump() const
Returns the const pointer to receiver&#39;s temporary jump.
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
const FloatMatrix & giveF() const
Returns the const pointer to receiver&#39;s deformation gradient vector.
const FloatArray & giveNormal() const
Returns const reference to normal vector.
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
Class representing solution step.
Definition: timestep.h:80
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in Reduced form.
Definition: material.C:142

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