OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
intelline2intpen.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 "intelline2intpen.h"
36 
37 #include "dofmanager.h"
38 #include "fei2dlinelin.h"
39 #include "fei2dlinequad.h"
40 #include "gausspoint.h"
41 #include "gaussintegrationrule.h"
42 
43 
44 namespace oofem {
45 REGISTER_Element(IntElLine2IntPen);
46 
48  numberOfDofMans = 6;
49 
51 }
52 
54 
55 }
56 
59 {
60  this->axisymmode = false;
63 
64 
65  // Check if node numbering is ok
66 // printf("this->dofManArray: "); this->dofManArray.printYourself();
67 // DofManager *node1 = this->giveDofManager(1);
68 // const FloatArray &x1 = *(node1->giveCoordinates());
69 // printf("x1: "); x1.printYourself();
70 
71 // DofManager *node2 = this->giveDofManager(2);
72 // const FloatArray &x2 = *(node2->giveCoordinates());
73 // printf("x2: "); x2.printYourself();
74 
75 // DofManager *node3 = this->giveDofManager(3);
76 // const FloatArray &x3 = *(node3->giveCoordinates());
77 // printf("x3: "); x3.printYourself();
78 
79 // DofManager *node4 = this->giveDofManager(4);
80 // const FloatArray &x4 = *(node4->giveCoordinates());
81 // printf("x4: "); x4.printYourself();
82 
83 // IntArray nodeInd = {dofManArray.at(2), dofManArray.at(3), dofManArray.at(4)};
84 
85 // double L2 = x1.distance_square(x2);
86 // printf("L2: %e\n", L2);
87 
88 // double L3 = x1.distance_square(x3);
89 // printf("L3: %e\n", L3);
90 
91 // double L4 = x1.distance_square(x4);
92 // printf("L4: %e\n", L4);
93 //
94 // FloatArray dist2 = {L2, L3, L4};
95 //
96 // double minDist2 = 1.0e20, maxDist2 =-1.0;
97 // int minDistInd = -1, maxDistInd = -1;
98 
99 // for(int i = 1; i <= 3; i++ ) {
100 // if(dist2.at(i) <= minDist2) {
101 // minDist2 = dist2.at(i);
102 // minDistInd = i;
103 // }
104 //
105 // if(dist2.at(i) >= maxDist2) {
106 // maxDist2 = dist2.at(i);
107 // maxDistInd = i;
108 // }
109 // }
110 //
111 // printf("minDistInd: %d\n", minDistInd);
112 // printf("maxDistInd: %d\n", maxDistInd);
113 
114 
115 // int midDistInd = -1;
116 // for(int i = 1; i <= 3; i++ ) {
117 // if(i != minDistInd && i != maxDistInd) {
118 // midDistInd = i;
119 // }
120 // }
121 // printf("midDistInd: %d\n", midDistInd);
122 
123 
124 // if(L2 < L3) {
125 // printf("Renumbering element %d\n.\n", this->giveNumber());
126 // std::swap(dofManArray.at(2), dofManArray.at(3));
127 // }
128 
129 
130  return result;
131 }
132 
133 
134 void
136 {
137 // printf("Entering IntElLine2IntPen :: computeCovarBaseVectorAt\n");
138 
139  // Since we are averaging over the whole element, always evaluate the base vectors at xi = 0.
140 
141  FloatArray xi_0 = {0.0};
142 
143  FloatMatrix dNdxi;
145 // interp->evaldNdxi( dNdxi, ip->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
146  interp->evaldNdxi( dNdxi, xi_0, FEIElementGeometryWrapper(this) );
147 
148  G.resize(2);
149  G.zero();
150  int numNodes = this->giveNumberOfNodes();
151  for ( int i = 1; i <= dNdxi.giveNumberOfRows(); i++ ) {
152  double X1_i = 0.5 * ( this->giveNode(i)->giveCoordinate(1) + this->giveNode(i + numNodes / 2)->giveCoordinate(1) ); // (mean) point on the fictious mid surface
153  double X2_i = 0.5 * ( this->giveNode(i)->giveCoordinate(2) + this->giveNode(i + numNodes / 2)->giveCoordinate(2) );
154  G.at(1) += dNdxi.at(i, 1) * X1_i;
155  G.at(2) += dNdxi.at(i, 1) * X2_i;
156  }
157 }
158 
159 void
161 {
162  // Computes the stiffness matrix of the receiver K_cohesive = int_A ( N^t * dT/dj * N ) dA
163  FloatMatrix N, D, DN;
164  bool matStiffSymmFlag = this->giveCrossSection()->isCharacteristicMtrxSymmetric(rMode);
165  answer.clear();
166 
167  FloatMatrix rotationMatGtoL;
168  FloatArray u;
169 
170  // First loop over GP: compute projection of test function and traction.
171  // The setting is as follows: we have an interface with quadratic interpolation and we
172  // wish to project onto the space of constant functions on each element.
173 
174  // Projecting the basis functions gives a constant for each basis function.
175  FloatMatrix proj_N;
176  proj_N.clear();
177 
178  FloatMatrix proj_DN;
179  proj_DN.clear();
180 
181  double area = 0.;
182 
183  for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {
184 
185  this->computeNmatrixAt(ip, N);
186 
187  if ( this->nlGeometry == 0 ) {
188  this->giveStiffnessMatrix_Eng(D, rMode, ip, tStep);
189  } else if ( this->nlGeometry == 1 ) {
190  this->giveStiffnessMatrix_dTdj(D, rMode, ip, tStep);
191  } else {
192  OOFEM_ERROR("nlgeometry must be 0 or 1!")
193  }
194 
195  this->computeTransformationMatrixAt(ip, rotationMatGtoL);
196  D.rotatedWith(rotationMatGtoL, 'n'); // transform stiffness to global coord system
197 
198  this->computeNmatrixAt(ip, N);
199  DN.beProductOf(D, N);
200 
201 
202  double dA = this->computeAreaAround(ip);
203  area += dA;
204 
205  proj_N.add(dA, N);
206  proj_DN.add(dA, DN);
207  }
208 
209 // printf("area: %e\n", area);
210  proj_N.times(1./area);
211  proj_DN.times(1./area);
212 
213 // printf("proj_N: "); proj_N.printYourself();
214 
215  for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {
216 
217  if ( this->nlGeometry == 0 ) {
218  this->giveStiffnessMatrix_Eng(D, rMode, ip, tStep);
219  } else if ( this->nlGeometry == 1 ) {
220  this->giveStiffnessMatrix_dTdj(D, rMode, ip, tStep);
221  } else {
222  OOFEM_ERROR("nlgeometry must be 0 or 1!")
223  }
224 
225  this->computeTransformationMatrixAt(ip, rotationMatGtoL);
226  D.rotatedWith(rotationMatGtoL, 'n'); // transform stiffness to global coord system
227 
228  DN.beProductOf(D, proj_N);
229 
230 
231  double dA = this->computeAreaAround(ip);
232  if ( matStiffSymmFlag ) {
233  answer.plusProductSymmUpper(proj_N, DN, dA);
234  } else {
235  answer.plusProductUnsym(proj_N, DN, dA);
236  }
237  }
238 
239 
240  if ( matStiffSymmFlag ) {
241  answer.symmetrized();
242  }
243 }
244 
245 
246 void
248  TimeStep *tStep, int useUpdatedGpRecord)
249 {
250  // Computes internal forces
251  // For this element we use an "interior penalty" formulation, where
252  // the cohesive zone contribution is weakened, i.e. the traction and
253  // test function for the cohesive zone are projected onto a reduced
254  // space. The purpose of the projection is to improve the stability
255  // properties of the formulation, thereby avoiding traction oscilations.
256 
257  FloatMatrix N;
258  FloatArray u, traction, jump;
259 
260  this->computeVectorOf(VM_Total, tStep, u);
261  // subtract initial displacements, if defined
262  if ( initialDisplacements.giveSize() ) {
264  }
265 
266  // zero answer will resize accordingly when adding first contribution
267  answer.clear();
268 
269 
270 
271  // First loop over GP: compute projection of test function and traction.
272  // The setting is as follows: we have an interface with quadratic interpolation and we
273  // wish to project onto the space of constant functions on each element.
274 
275  // The projection of t becomes a constant
276 // FloatArray proj_t;
277 // proj_t.clear();
278 
279 
280  FloatArray proj_jump;
281  proj_jump.clear();
282 
283  // Projecting the basis functions gives a constant for each basis function.
284  FloatMatrix proj_N;
285  proj_N.clear();
286 
287  double area = 0.;
288 
289  for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {
290 
291  this->computeNmatrixAt(ip, N);
292  jump.beProductOf(N, u);
293 // this->computeTraction(traction, ip, jump, tStep);
294 
295  double dA = this->computeAreaAround(ip);
296  area += dA;
297 
298  proj_jump.add(dA, jump);
299  proj_N.add(dA, N);
300  }
301 
302 // printf("area: %e\n", area);
303  proj_jump.times(1./area);
304  proj_N.times(1./area);
305 
306 // printf("proj_N: "); proj_N.printYourself();
307 
308 
309  // Second loop over GP: assemble contribution to internal forces as we are used to.
310  for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {
311 // this->computeNmatrixAt(ip, N);
312 // jump.beProductOf(N, u);
313  this->computeTraction(traction, ip, proj_jump, tStep);
314 
315  // compute internal cohesive forces as f = N^T*traction dA
316  double dA = this->computeAreaAround(ip);
317  answer.plusProduct(proj_N, traction, dA);
318  }
319 
320 }
321 
322 
323 
324 } /* namespace oofem */
CrossSection * giveCrossSection()
Definition: element.C:495
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
This class implements a two dimensional interface element and is simply an extension of IntElLine1 to...
Definition: intelline2.h:52
virtual void evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
Definition: feinterpol.h:193
Class and object Domain.
Definition: domain.h:115
void computeVectorOf(ValueModeType u, TimeStep *tStep, FloatArray &answer)
Returns local vector of unknowns.
Definition: element.C:86
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
Access method for default integration rule.
Definition: element.h:822
IntElLine2IntPen(int n, Domain *d)
bool axisymmode
Flag controlling axisymmetric mode (integration over unit circumferential angle)
Definition: intelline1.h:60
virtual void computeTraction(FloatArray &traction, IntegrationPoint *ip, const FloatArray &jump, TimeStep *tStep)
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
virtual void computeTransformationMatrixAt(GaussPoint *gp, FloatMatrix &answer)
Definition: intelline1.C:197
virtual void computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
Computes the stiffness/tangent matrix of receiver.
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual double giveCoordinate(int i)
Definition: node.C:82
virtual double computeAreaAround(GaussPoint *gp)
Definition: intelline1.C:121
MatResponseMode
Describes the character of characteristic material matrix.
FloatArray initialDisplacements
Initial displacement vector, describes the initial nodal displacements when element has been casted...
virtual int giveNumberOfNodes() const
Returns number of nodes of receiver.
Definition: element.h:662
Class representing a general abstraction for finite element interpolation class.
Definition: feinterpol.h:132
void plusProductUnsym(const FloatMatrix &a, const FloatMatrix &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:780
#define OOFEM_ERROR(...)
Definition: error.h:61
REGISTER_Element(LSpace)
virtual void giveStiffnessMatrix_Eng(FloatMatrix &answer, MatResponseMode rMode, IntegrationPoint *ip, TimeStep *tStep)
Definition: intelline1.h:89
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
void rotatedWith(const FloatMatrix &r, char mode= 'n')
Returns the receiver &#39;a&#39; transformed using give transformation matrix r.
Definition: floatmatrix.C:1557
Wrapper around element definition to provide FEICellGeometry interface.
Definition: feinterpol.h:95
virtual bool isCharacteristicMtrxSymmetric(MatResponseMode rMode)
Check for symmetry of stiffness matrix.
Definition: crosssection.h:172
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
#define N(p, q)
Definition: mdm.C:367
void plusProductSymmUpper(const FloatMatrix &a, const FloatMatrix &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:698
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
#define _IFT_IntElLine1_axisymmode
Definition: intelline1.h:41
virtual void giveStiffnessMatrix_dTdj(FloatMatrix &answer, MatResponseMode rMode, IntegrationPoint *ip, TimeStep *tStep)
int numberOfGaussPoints
Number of integration points as specified by nip.
Definition: element.h:188
virtual void computeCovarBaseVectorAt(GaussPoint *gp, FloatArray &G)
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual FEInterpolation * giveInterpolation() const
Definition: intelline2.C:114
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
Class representing the general Input Record.
Definition: inputrecord.h:101
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
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
virtual void giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord=0)
Returns equivalent nodal forces vectors.
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
static FEI2dLineQuad interp
Definition: intelline2.h:55
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.
void clear()
Sets size of receiver to be an empty matrix. It will have zero rows and zero columns size...
Definition: floatmatrix.h:516
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
void symmetrized()
Initializes the lower half of the receiver according to the upper half.
Definition: floatmatrix.C:1576
int nlGeometry
Flag indicating if geometrical nonlinearities apply.
Class representing integration point in finite element program.
Definition: gausspoint.h:93
Class representing solution step.
Definition: timestep.h:80
int numberOfDofMans
Number of dofmanagers.
Definition: element.h:149
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
virtual void computeNmatrixAt(GaussPoint *gp, FloatMatrix &answer)
Computes modified interpolation matrix (N) for the element which multiplied with the unknowns vector ...
Definition: intelline2.C:68
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:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011