OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
contactelement.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 "contact/contactelement.h"
36 #include "dofmanager.h"
37 #include "floatarray.h"
38 #include "floatmatrix.h"
39 #include "valuemodetype.h"
40 #include "dofiditem.h"
41 #include "timestep.h"
42 #include "dof.h"
43 #include "masterdof.h"
44 #include "unknownnumberingscheme.h"
45 #include "domain.h"
46 #include "gaussintegrationrule.h"
47 #include "error.h"
48 
49 namespace oofem {
50 
52 {
53  this->dofIdArray.clear();
54  this->integrationRule = NULL;
55 };
56 
57 
58 #if 0
60 {
61  this->masterNode = master;
62  this->slaveNode = slave;
63  this->area = 1.0; // should be optional parameter
64  this->epsN = 1.0e6; // penalty - should be given by 'contactmaterial'
65 };
66 
67 int
69 {
70  // compute normal as direction vector from master node to slave node
71  FloatArray xs, xm, normal;
72  xs = *this->slaveNode->giveCoordinates();
73  xm = *this->masterNode->giveCoordinates();
74 
75  normal = xs-xm;
76  double norm = normal.computeNorm();
77  if ( norm < 1.0e-8 ) {
78  OOFEM_ERROR("Couldn't compute normal between master node (num %d) and slave node (num %d), nodes are too close to each other.",
79  masterNode->giveGlobalNumber(), slaveNode->giveGlobalNumber() )
80  } else {
81  this->normal = normal*(1.0/norm);
82  }
83  return 1;
84 }
85 
86 
87 void
89 {
90 
91  FloatArray xs, xm, uS, uM;
92  xs = *this->slaveNode->giveCoordinates();
93  xm = *this->masterNode->giveCoordinates();
94  this->slaveNode->giveUnknownVector(uS, {D_u, D_v, D_w}, VM_Total, tStep, true);
95  this->masterNode->giveUnknownVector(uM, {D_u, D_v, D_w}, VM_Total, tStep, true);
96  xs.add(uS);
97  xm.add(uM);
98  FloatArray dx = xs-xm;
99 
100  FloatArray normal = this->giveNormal();
101  answer = {dx.dotProduct(normal), 0.0, 0.0};
102 
103  printf("normal gap = %e \n", answer.at(1));
104  if ( answer.at(1) < 0.0 ) {
105  //printf("normal gap = %e \n", answer.at(1));
106  this->inContact = true; // store in gp?
107  } else {
108  this->inContact = false;
109  }
110 
111 }
112 
113 
114 void
116 {
117 
118  // The normal is not updated for node2node which is for small deformations only
119  // C = {n -n}
120  FloatArray normal = this->giveNormal();
121  answer = { normal.at(1), normal.at(2), normal.at(3),
122  -normal.at(1), -normal.at(2), -normal.at(3) };
123 
124 }
125 
126 
127 
128 void
130 {
131  // should be replaced with a call to constitutive model
132  // gap should be in a local system
133  if ( gap.at(1) < 0.0 ) {
134  t = this->epsN * gap;
135  } else {
136  t = {0.0, 0.0, 0.0};
137  }
138 
139 }
140 
141 
142 
143 
144 
145 void
147  const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms)
148 {
149  answer.clear();
150  FloatArray gap, C;
151 
152  this->computeGap(gap, tStep);
153  if ( gap.at(1) < 0.0 ) {
155  FloatArray t;
156  this->computeContactTractionAt(gp, t ,gap, tStep);
157 
158  this->computeCmatrixAt(gp, C, tStep);
159 
160  // compute load vector
161  // fc = C^T * traction * A, Area - should be optional par
162  answer = t.at(1) * this->area * C;
163 
164  }
165 
166 }
167 
168 
169 void
170 Node2NodeContact :: computeContactTangent(FloatMatrix &answer, CharType type, TimeStep *tStep)
171 {
172  // Need to set up an integration rule
173  GaussPoint *gp = NULL;
174  FloatArray gap;
175 
176  this->computeGap(gap, tStep);
177 
178  FloatArray C;
179  this->computeCmatrixAt(gp, C, tStep);
180  answer.beDyadicProductOf(C,C);
181  // this is the interface stiffness and should be obtained from that model
182  answer.times( this->epsN * this->area );
183  answer.negated();
184 
185  if( gap.at(1) > 0.0 ) {
186  answer.zero();
187  }
188 }
189 
190 
191 
192 
193 
194 
195 void
197 {
198  // should return location array for the master and the slave node
199  // TODO this whole thing should be rewritten
200 
201 
202  answer.resize(6);
203  answer.zero();
204  //TODO should use a proper unknownnumberingscheme
205  IntArray dofIdArray = {D_u, D_v, D_w};
206 
207  // master node
208  for ( int i = 1; i <= dofIdArray.giveSize(); i++ ) {
209  if ( this->masterNode->hasDofID( (DofIDItem)dofIdArray.at(i) ) ) { // add corresponding number
210  Dof *dof= this->masterNode->giveDofWithID( (DofIDItem)dofIdArray.at(i) );
211  answer.at(i) = s.giveDofEquationNumber( dof );
212  }
213  }
214 
215  // slave node
216  for ( int i = 1; i <= dofIdArray.giveSize(); i++ ) {
217  if ( this->slaveNode->hasDofID( (DofIDItem)dofIdArray.at(i) ) ) { // add corresponding number
218  Dof *dof= this->slaveNode->giveDofWithID( (DofIDItem)dofIdArray.at(i) );
219  answer.at(3 + i) = s.giveDofEquationNumber( dof );
220  }
221  }
222 
223 }
224 
225 
226 void
228 {
229  // Sets up the integration rule array which contains all the necessary integration points
230  if ( this->integrationRule == NULL ) {
231  //TODO sets a null pointer for the element in the iRule
232  this->integrationRule = new GaussIntegrationRule(1, NULL) ;
233  this->integrationRule->SetUpPoint(_Unknown);
234  }
235 
236 
237 }
238 
239 
240 
241 
242 
243 
244 
245 
246 
247 
248 
249 
250 
251 // node 2 node Lagrange
252 
253 
255 {
256  this->masterNode = master;
257  this->slaveNode = slave;
258  this->area = 1.0;
259 };
260 
261 
262 
263 void
265 {
266 
268 
269  // Add one lagrange dof
270  if ( this->masterNode->hasDofID( (DofIDItem)this->giveDofIdArray().at(1) ) ) {
271  Dof *dof= this->masterNode->giveDofWithID( (DofIDItem)this->giveDofIdArray().at(1) );
272  answer.followedBy( s.giveDofEquationNumber(dof) );
273  }
274 
275 }
276 
277 
278 
279 void
281  const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms)
282 {
283 
284  //Loop through all the master objects and let them do their thing
285  FloatArray gap, C, Fc;
286  this->computeGap(gap, tStep);
287 
289  FloatArray t;
290  this->computeContactTractionAt(gp, t ,gap, tStep);
291  this->computeCmatrixAt(gp, C, tStep);
292 
293  // compute load vector
294  // for Lagrange: fc = traction * C^T * A (traction = lambda)
295  FloatArray temp = t.at(1) *this->area * C;
296 
297  answer.resize( C.giveSize() + 1);
298  answer.zero();
299  if( gap.at(1) < 0.0 ) {
300  answer.addSubVector(temp,1);
301  answer.at( C.giveSize() + 1 ) = -gap.at(1);
302  }
303 }
304 
305 
306 
307 void
309 {
310  answer.resize(7,7);
311  answer.zero();
312 
313  FloatArray gap;
314  this->computeGap(gap, tStep);
315 
316  if( gap.at(1) < 0.0 ) {
317 
319 
320  FloatArray C;
321  this->computeCmatrixAt(gp, C, tStep);
322  int sz = C.giveSize();
323  C.times(this->area);
324 
325  answer.addSubVectorCol(C, 1, sz + 1);
326  answer.addSubVectorRow(C, sz + 1, 1);
327  }
328 
329  //TODO need to add a small number for the solver
330  for ( int i = 1; i <= 7; i++ ) {
331  answer.at(i,i) += 1.0e-8;
332  }
333 
334 }
335 
336 
337 
338 void
340 {
341  // should be replaced with a call to constitutive model
342  // gap should be in a local system
343  if ( gap.at(1) < 0.0 ) {
344 
345  Dof *dof = masterNode->giveDofWithID( this->giveDofIdArray().at(1) );
346  double lambda = dof->giveUnknown(VM_Total, tStep);
347  t = {lambda, 0.0, 0.0};
348  printf("lambda %e \n\n", lambda);
349  } else {
350  t = {0.0, 0.0, 0.0};
351  }
352 
353 }
354 
355 
356 void
358 {
359  answer = {this->masterNode->giveNumber()};
360 }
361 
362 #endif
363 
364 }
365 
366 
367 
368 
369 
370 
371 
372 
373 
374 
375 
376 
virtual void giveLocationArray(IntArray &answer, const UnknownNumberingScheme &s)
Definition: celnode2node.C:163
int SetUpPoint(MaterialMode mode)
Trivial implementation, only creates a single point.
Class and object Domain.
Definition: domain.h:115
virtual void computeContactTractionAt(GaussPoint *gp, FloatArray &t, FloatArray &gap, TimeStep *tStep)
Definition: celnode2node.C:281
ContactElement()
Constructor.
void zero()
Sets all component to zero.
Definition: intarray.C:52
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
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
virtual double giveUnknown(ValueModeType mode, TimeStep *tStep)=0
The key method of class Dof.
virtual IntArray & giveDofIdArray()
virtual int instanciateYourself(DataReader &dr)
Definition: celnode2node.C:53
Base class for dof managers.
Definition: dofmanager.h:113
virtual void computeContactTangent(FloatMatrix &answer, TimeStep *tStep)
Definition: celnode2node.C:141
void negated()
Changes sign of receiver values.
Definition: floatmatrix.C:1605
Class representing the abstraction for input data source.
Definition: datareader.h:50
virtual void computeCmatrixAt(GaussPoint *gp, FloatArray &answer, TimeStep *TimeStep)
Definition: celnode2node.C:94
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
virtual void giveDofManagersToAppendTo(IntArray &answer)
Definition: celnode2node.C:297
virtual int giveDofEquationNumber(Dof *dof) const =0
Returns the equation number for corresponding DOF.
Node2NodeContactL(DofManager *master, DofManager *slave)
Constructor.
Definition: celnode2node.C:204
double dotProduct(const FloatArray &x) const
Computes the dot product (or inner product) of receiver and argument.
Definition: floatarray.C:463
#define OOFEM_ERROR(...)
Definition: error.h:61
void clear()
Clears the array (zero size).
Definition: intarray.h:177
DofIDItem
Type representing particular dof type.
Definition: dofiditem.h:86
Node2NodeContact(DofManager *master, DofManager *slave)
Constructor.
Definition: celnode2node.C:43
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
Abstract base class allowing to control the way, how equations are assigned to individual DOFs...
void addSubVector(const FloatArray &src, int si)
Adds the given vector as sub-vector to receiver.
Definition: floatarray.C:399
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
virtual void giveLocationArray(IntArray &answer, const UnknownNumberingScheme &s)
Definition: celnode2node.C:213
virtual void computeGap(FloatArray &answer, TimeStep *tStep)
Definition: celnode2node.C:73
virtual void setupIntegrationPoints()
Definition: celnode2node.C:192
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
GaussPoint * getIntegrationPoint(int n)
Access particular integration point of receiver.
void addSubVectorCol(const FloatArray &src, int sr, int sc)
Adds given vector to receiver starting at given position.
Definition: floatmatrix.C:627
double norm(const FloatArray &x)
Definition: floatarray.C:985
CharType
Definition: chartype.h:87
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual void computeContactTangent(FloatMatrix &answer, TimeStep *tStep)
Definition: celnode2node.C:252
void followedBy(const IntArray &b, int allocChunk=0)
Appends array b at the end of receiver.
Definition: intarray.C:145
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
void beDyadicProductOf(const FloatArray &vec1, const FloatArray &vec2)
Assigns to the receiver the dyadic product .
Definition: floatmatrix.C:492
IntegrationRule * integrationRule
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
int giveSize() const
Definition: intarray.h:203
virtual void computeContactForces(FloatArray &answer, TimeStep *tStep, ValueModeType mode, const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms)
Definition: celnode2node.C:226
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
the oofem namespace is to define a context or scope in which all oofem names are defined.
Abstract class Dof represents Degree Of Freedom in finite element mesh.
Definition: dof.h:93
virtual void computeContactForces(FloatArray &answer, TimeStep *tStep, ValueModeType mode, const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms)
Definition: celnode2node.C:118
void addSubVectorRow(const FloatArray &src, int sr, int sc)
Adds given vector to receiver starting at given position.
Definition: floatmatrix.C:607
Class representing integration point in finite element program.
Definition: gausspoint.h:93
virtual void computeContactTractionAt(GaussPoint *gp, FloatArray &t, FloatArray &gap, TimeStep *tStep)
Definition: celnode2node.C:105
Class representing solution step.
Definition: timestep.h:80
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
Class representing Gaussian-quadrature integration rule.
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:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011