1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbutil/gauss.cc,v 1.23 2017/01/12 14:44:04 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31
32 /* Punti di Gauss */
33
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35
36 #include "gauss.h"
37
38
39 /* Dati dei punti */
40 const doublereal dGaussPnt[] = {
41 /* 1 point, 2nd order */ 0.000000000000000,
42 /* 2 point, 4th order */ -0.577350269189626, 0.577350269189626,
43 /* 3 point, 6th order */ -0.774596669241483, 0.000000000000000, 0.774596669241483,
44 /* 4 point, 8th order */ -0.861136311594053, -0.339981043584856, 0.339981043584856, 0.861136311594053,
45 /* 5 point, 10th order */ -0.906179845938664, -0.538469310105683, 0.000000000000000, 0.538469310105683, 0.906179845938664,
46 /* 6 point, 12th order */ -0.932469514203152, -0.661209386466264, -0.238619186083197, 0.238619186083197, 0.661209386466264, 0.932469514203152,
47 /* 7 point, 14th order */ -0.949107912342758, -0.741531185599394, -0.405845151377397, 0.000000000000000, 0.405845151377397, 0.741531185599394, 0.949107912342758,
48 /* 8 point, 16th order */ -0.960289856497536, -0.796666477413627, -0.525532409916329, -0.183434642495650, 0.183434642495650, 0.525532409916329, 0.796666477413627, 0.960289856497536,
49 /* 9 point, 18th order */ -0.968160239507626, -0.836031107326636, -0.613371432700590, -0.324253423403809, 0.000000000000000, 0.324253423403809, 0.613371432700591, 0.836031107326636, 0.968160239507626,
50 /* 10 point, 20th order */ -0.973906528517171, -0.865063366688985, -0.679409568299025, -0.433395394129247, -0.148874338981631, 0.148874338981631, 0.433395394129247, 0.679409568299024, 0.865063366688984, 0.973906528517172
51 /* 11 point, 22nd order */
52 };
53
54 const doublereal dGaussWght[] = {
55 /* 1 point, 2nd order */ 2.000000000000000,
56 /* 2 point, 4nd order */ 1.000000000000000, 1.000000000000000,
57 /* 3 point, 6nd order */ 0.555555555555556, 0.888888888888889, 0.555555555555556,
58 /* 4 point, 8nd order */ 0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454,
59 /* 5 point, 10nd order */ 0.236926885056189, 0.478628670499367, 0.568888888888889, 0.478628670499367, 0.236926885056189,
60 /* 6 point, 12th order */ 0.171324492379171, 0.360761573048139, 0.467913934572690, 0.467913934572690, 0.360761573048139, 0.171324492379171,
61 /* 7 point, 14th order */ 0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277, 0.129484966168870,
62 /* 8 point, 16th order */ 0.101228536290376, 0.222381034453375, 0.313706645877887, 0.362683783378362, 0.362683783378362, 0.313706645877887, 0.222381034453375, 0.101228536290376,
63 /* 9 point, 18th order */ 0.081274388361575, 0.180648160694858, 0.260610696402935, 0.312347077040002, 0.330239355001259, 0.312347077040002, 0.260610696402936, 0.180648160694857, 0.081274388361575,
64 /* 10 point, 20th order */ 0.066671344308688, 0.149451349150581, 0.219086362515981, 0.269266719309996, 0.295524224714752, 0.295524224714752, 0.269266719309996, 0.219086362515982, 0.149451349150581, 0.066671344308688
65 /* 11 point, 22nd order */
66 };
67
68 const doublereal dTrapezoidPnt[] = {
69 0.,
70 -1., 1.,
71 -1., 0., 1.,
72 -1., -.3333333333333333, .3333333333333333, 1.,
73 -1., -.5, 0., .5, 1.
74 };
75
76 const doublereal dTrapezoidWght[] = {
77 2.,
78 1., 1.,
79 .3333333333333333, 1.333333333333333, .3333333333333333,
80 .25, .75, .75, .25,
81 .1555555555555556, .7111111111111111, .2666666666666667,
82 .7111111111111111, .1555555555555556
83 };
84
85 /* GaussData - begin */
86
GaussData(integer iN)87 GaussData::GaussData(integer iN)
88 : NumIntData((iN < 1 ? 1 : (iN > 10 ? 10 : iN))), pdPnt(NULL), pdWght(NULL) {
89 ASSERT(iNum >= 1 && iNum <= 10);
90 integer i = iNum*(iNum-1)/2-1;
91 (doublereal*&)pdPnt = (doublereal*)dGaussPnt+i;
92 (doublereal*&)pdWght = (doublereal*)dGaussWght+i;
93 }
94
95
dGetPnt(integer i) const96 doublereal GaussData::dGetPnt(integer i) const
97 {
98 ASSERT(i > 0 && i <= iNum);
99 return pdPnt[i];
100 }
101
102
dGetWght(integer i) const103 doublereal GaussData::dGetWght(integer i) const
104 {
105 ASSERT(i > 0 && i <= iNum);
106 return pdWght[i];
107 }
108
109
Get(integer i) const110 PntWght GaussData::Get(integer i) const
111 {
112 ASSERT(i > 0 && i <= iNum);
113 return PntWght(pdPnt[i], pdWght[i]);
114 }
115
116
pdGetPnt(void) const117 const doublereal* GaussData::pdGetPnt(void) const
118 {
119 return pdPnt+1;
120 }
121
122
pdGetWght(void) const123 const doublereal* GaussData::pdGetWght(void) const
124 {
125 return pdWght+1;
126 }
127
128 /* GaussData - end */
129
130
131 /* TrapezoidData - begin */
132
TrapezoidData(integer iN)133 TrapezoidData::TrapezoidData(integer iN)
134 : NumIntData((iN < 1 ? 1 : (iN > 5 ? 5 : iN))), pdPnt(NULL), pdWght(NULL) {
135 integer i = iNum*(iNum-1)/2-1;
136 (doublereal*&)pdPnt = (doublereal*)dTrapezoidPnt+i;
137 (doublereal*&)pdWght = (doublereal*)dTrapezoidWght+i;
138 }
139
140
dGetPnt(integer i) const141 doublereal TrapezoidData::dGetPnt(integer i) const
142 {
143 ASSERT(i > 0 && i <= iNum);
144 return pdPnt[i];
145 }
146
147
dGetWght(integer i) const148 doublereal TrapezoidData::dGetWght(integer i) const
149 {
150 ASSERT(i > 0 && i <= iNum);
151 return pdWght[i];
152 }
153
154
Get(integer i) const155 PntWght TrapezoidData::Get(integer i) const
156 {
157 ASSERT(i > 0 && i <= iNum);
158 return PntWght(pdPnt[i], pdWght[i]);
159 }
160
161
pdGetPnt(void) const162 const doublereal* TrapezoidData::pdGetPnt(void) const
163 {
164 return pdPnt+1;
165 }
166
167
pdGetWght(void) const168 const doublereal* TrapezoidData::pdGetWght(void) const
169 {
170 return pdWght+1;
171 }
172
173 /* TrapezoidData - end */
174
175
176 /* GaussDataIterator - begin */
177
GaussDataIterator(integer iN)178 GaussDataIterator::GaussDataIterator(integer iN)
179 : GaussData(iN), iCurr(1)
180 {
181 NO_OP;
182 }
183
184
dGetFirst(integer i) const185 doublereal GaussDataIterator::dGetFirst(integer i) const
186 {
187 ASSERT(i == 0 || i == 1);
188
189 (integer&)iCurr = 1;
190 if(i == 0) {
191 return pdPnt[1];
192 }
193 /* if(i == 1) */
194 return pdWght[1];
195 }
196
197
GetFirst(void) const198 PntWght GaussDataIterator::GetFirst(void) const
199 {
200 (integer&)iCurr = 1;
201 return PntWght(pdPnt[1], pdWght[1]);
202 }
203
204
fGetNext(doublereal & d,integer i) const205 flag GaussDataIterator::fGetNext(doublereal& d, integer i) const
206 {
207 ASSERT(i == 0 || i == 1);
208
209 (integer&)iCurr += 1;
210 if(iCurr > iNum) {
211 return flag(0);
212 }
213
214 if(i == 0) {
215 d = pdPnt[iCurr];
216 } else /* if(i == 1) */ {
217 d = pdWght[iCurr];
218 }
219
220 return flag(1);
221 }
222
223
dGetCurrPnt(void) const224 doublereal GaussDataIterator::dGetCurrPnt(void) const
225 {
226 ASSERT(iCurr >= 1 && iCurr <= iNum);
227
228 return pdPnt[iCurr];
229 }
230
231
dGetCurrWght(void) const232 doublereal GaussDataIterator::dGetCurrWght(void) const
233 {
234 ASSERT(iCurr >= 1 && iCurr <= iNum);
235
236 return pdWght[iCurr];
237 }
238
239
fGetNext(PntWght & PW) const240 flag GaussDataIterator::fGetNext(PntWght& PW) const
241 {
242 (integer&)iCurr += 1;
243 if(iCurr > iNum) {
244 return flag(0);
245 }
246
247 PW = PntWght(pdPnt[iCurr], pdWght[iCurr]);
248
249 return flag(1);
250 }
251
252 /* GaussdataIterator - end */
253
254
255 /* NumIntIterator - begin */
256
NumIntIterator(NumIntData & d)257 NumIntIterator::NumIntIterator(NumIntData& d)
258 : iCurr(1), data(d)
259 {
260 NO_OP;
261 }
262
263
~NumIntIterator(void)264 NumIntIterator::~NumIntIterator(void)
265 {
266 NO_OP;
267 }
268
269
dGetFirst(integer i) const270 doublereal NumIntIterator::dGetFirst(integer i) const
271 {
272 ASSERT(i == 0 || i == 1);
273
274 (integer&)iCurr = 1;
275 if(i == 0) {
276 return data.dGetPnt(1);
277 }
278 /* else if(i == 1) */
279 return data.dGetWght(1);
280 }
281
282
GetFirst(void) const283 PntWght NumIntIterator::GetFirst(void) const
284 {
285 (integer&)iCurr = 1;
286 return PntWght(data.Get(1));
287 }
288
289
fGetNext(doublereal & d,integer i) const290 flag NumIntIterator::fGetNext(doublereal& d, integer i) const
291 {
292 ASSERT(i == 0 || i == 1);
293
294 (integer&)iCurr += 1;
295 if(iCurr > data.iGetNum()) {
296 return flag(0);
297 }
298
299 if(i == 0) {
300 d = data.dGetPnt(iCurr);
301 } else /* if(i == 1) */ {
302 d = data.dGetWght(iCurr);
303 }
304
305 return flag(1);
306 }
307
308
dGetCurrPnt(void) const309 doublereal NumIntIterator::dGetCurrPnt(void) const
310 {
311 ASSERT(iCurr >= 1 && iCurr <= data.iGetNum());
312
313 return data.dGetPnt(iCurr);
314 }
315
316
dGetCurrWght(void) const317 doublereal NumIntIterator::dGetCurrWght(void) const
318 {
319 ASSERT(iCurr >= 1 && iCurr <= data.iGetNum());
320
321 return data.dGetWght(iCurr);
322 }
323
324
fGetNext(PntWght & PW) const325 flag NumIntIterator::fGetNext(PntWght& PW) const
326 {
327 (integer&)iCurr += 1;
328 if(iCurr > data.iGetNum()) {
329 return flag(0);
330 }
331
332 PW = PntWght(data.Get(iCurr));
333
334 return flag(1);
335 }
336
337 /* NumIntIterator - end */
338
339
340