OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
lobattoir.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
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 "lobattoir.h"
36 #include "gausspoint.h"
37 #include "floatarray.h"
38 #include "mathfem.h"
39
40 namespace oofem {
42  int startIndx, int endIndx, bool dynamic) :
43  IntegrationRule(n, e, startIndx, endIndx, dynamic) { }
44
46  IntegrationRule(n, e) { }
47
49 { }
50
51
52 int
54 {
55  FloatArray coords_xi, weights;
56  this->giveLineCoordsAndWeights(nPoints, coords_xi, weights);
57  this->gaussPoints.resize( nPoints );
58
59  for ( int i = 1; i <= nPoints; i++ ) {
60  this->gaussPoints [ i - 1 ] = new GaussPoint(this, i, {coords_xi.at(i)}, weights.at ( i ), mode);
61  }
62
63  this->intdomain = _Line;
64  return this->giveNumberOfIntegrationPoints();
65 }
66
67
68 int
70 //GaussIntegrationRule :: SetUpPointsOnSquare(int nPoints_xi1, int nPoints_xi2, MaterialMode mode)
71 {
72  int nPoints_xi1 = ( int ) floor( sqrt( double ( nPoints ) ) );
73  int nPoints_xi2 = nPoints_xi1;
74  FloatArray coords_xi1, weights1, coords_xi2, weights2;
75  this->giveLineCoordsAndWeights(nPoints_xi1, coords_xi1, weights1);
76  this->giveLineCoordsAndWeights(nPoints_xi2, coords_xi2, weights2);
77  this->gaussPoints.resize( nPoints_xi1 * nPoints_xi2 );
78  int count = 0;
79  for ( int i = 1; i <= nPoints_xi1; i++ ) {
80  for ( int j = 1; j <= nPoints_xi2; j++ ) {
81  count++;
82  this->gaussPoints [ count - 1 ] = new GaussPoint(this, count, {coords_xi1.at(i), coords_xi2.at(j)},
83  weights1.at ( i ) *weights2.at ( j ), mode);
84  }
85  }
86
87  this->intdomain = _Square;
88  return this->giveNumberOfIntegrationPoints();
89 }
90
91
92 int
94 //GaussIntegrationRule :: SetUpPointsOnCube(int nPoints_xi1, int nPoints_xi2, int nPoints_xi3, MaterialMode mode)
95 {
96  int nPoints_xi1 = ( int ) floor(cbrt( double ( nPoints ) ) + 0.5);
97  int nPoints_xi2 = nPoints_xi1;
98  int nPoints_xi3 = nPoints_xi1;
99  FloatArray coords_xi1, weights1, coords_xi2, weights2, coords_xi3, weights3;
100  this->giveLineCoordsAndWeights(nPoints_xi1, coords_xi1, weights1);
101  this->giveLineCoordsAndWeights(nPoints_xi2, coords_xi2, weights2);
102  this->giveLineCoordsAndWeights(nPoints_xi3, coords_xi3, weights3);
103  this->gaussPoints.resize( nPoints_xi1 * nPoints_xi2 * nPoints_xi3 );
104  int count = 0;
105  for ( int i = 1; i <= nPoints_xi1; i++ ) {
106  for ( int j = 1; j <= nPoints_xi2; j++ ) {
107  for ( int k = 1; k <= nPoints_xi3; k++ ) {
108  count++;
109  this->gaussPoints [ count - 1 ] = new GaussPoint(this, count, {coords_xi1.at(i), coords_xi2.at(j), coords_xi3.at(k)},
110  weights1.at ( i ) *weights2.at ( j ) *weights3.at ( k ), mode);
111  }
112  }
113  }
114
115  this->intdomain = _Cube;
116  return this->giveNumberOfIntegrationPoints();
117 }
118
119
120 int
122 {
123  OOFEM_ERROR("unsupported number of IPs (%d)", nPoints);
124  return 0;
125 }
126
127
128
129 int
131  int approxOrder)
132 {
133  int requiredNIP;
134
135  switch ( dType ) {
136  case _Line:
137  if ( approxOrder <= 1 ) {
138  return 1;
139  }
140
141  requiredNIP = ( int ) ceil( ( ( double ) approxOrder + 3.0 ) / 2. );
142  if ( requiredNIP > 6 ) {
143  return -1;
144  }
145
146  return requiredNIP;
147
148  default:
149  OOFEM_ERROR("unknown integrationDomain");
150  }
151
152  return -1;
153 }
154
155
156 void
158 // Create arrays of coordinates and weights for Lobatto Integration Points of a line with 'nPoints' integrationpoints
159 {
160  switch ( nPoints ) {
161  case 1:
162  coords_xi = FloatArray{0.0};
163  weights = FloatArray{2.0};
164  break;
165
166  case 2:
167  coords_xi = {-1.0, 1.0};
168  weights = {1.0, 1.0};
169  break;
170
171  case 3:
172  coords_xi = {-1.0, 0.0, 1.0};
173  weights = {0.333333333333333, 1.333333333333333, 0.333333333333333};
174  break;
175
176  case 4:
177  coords_xi = {-1.0,
178  -0.447213595499958,
179  0.447213595499958,
180  1.0};
181  weights = { 0.166666666666667,
182  0.833333333333333,
183  0.833333333333333,
184  0.166666666666667};
185  break;
186
187  case 5:
188  coords_xi = {-1.0,
189  -0.654653670707977,
190  0.0,
191  0.654653670707977,
192  1.0};
193  weights = { 0.1,
194  0.544444444444444,
195  0.711111111111111,
196  0.544444444444444,
197  0.1};
198  break;
199
200  default:
201  coords_xi = {-1.0,
202  -0.765055323929465,
203  -0.285231516480645,
204  0.285231516480645,
205  0.765055323929465,
206  1.0};
207  weights = { 0.066666666666667,
208  0.378474956297847,
209  0.554858377035486,
210  0.554858377035486,
211  0.378474956297847,
212  0.066666666666667};
213  if ( nPoints > 6 ) {
214  OOFEM_LOG_WARNING("Unsupported number of IPs (%d) for LobattoIR, using 6 ips instead.", nPoints);
215  }
216  break;
217  }
218 }
219 } // end namespace oofem
virtual int getRequiredNumberOfIntegrationPoints(integrationDomain dType, int approxOrder)
Abstract service.
Definition: lobattoir.C:130
virtual int SetUpPointsOnTriangle(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on triangular (area coords) integration domain.
Definition: lobattoir.C:121
integrationDomain
Used by integrator class to supply integration points for proper domain to be integrated (Area...
virtual ~LobattoIntegrationRule()
Destructor.
Definition: lobattoir.C:48
integrationDomain intdomain
Integration domain.
virtual int SetUpPointsOnSquare(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
Definition: lobattoir.C:69
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual int SetUpPointsOnCube(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on unit cube integration domain.
Definition: lobattoir.C:93
Abstract base class for all finite elements.
Definition: element.h:145
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
virtual int SetUpPointsOnLine(int nPoints, MaterialMode mode)
Sets up receiver&#39;s integration points on unit line integration domain.
Definition: lobattoir.C:53
Abstract base class representing integration rule.
#define OOFEM_LOG_WARNING(...)
Definition: logger.h:123
#define OOFEM_ERROR(...)
Definition: error.h:61
static void giveLineCoordsAndWeights(int nPoints, FloatArray &coords_xi, FloatArray &weights)
Definition: lobattoir.C:157
Class representing vector of real numbers.
Definition: floatarray.h:82
double cbrt(double x)
Returns the cubic root of x.
Definition: mathfem.h:109
int giveNumberOfIntegrationPoints() const
Returns number of integration points of receiver.
the oofem namespace is to define a context or scope in which all oofem names are defined.
Class representing integration point in finite element program.
Definition: gausspoint.h:93
std::vector< GaussPoint * > gaussPoints
Array containing integration points.
LobattoIntegrationRule(int n, Element *e, int startIndx, int endIndx, bool dynamic)
Constructor.
Definition: lobattoir.C:41