CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
cubtri.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine cubtri (F, T, EPS, MCALLS, ANS, ERR, NCALLS, W, NW, IDATA, RDATA, IER)
 
real *8 function rnderr (X, A, Y, B)
 
subroutine cubrul (F, VEC, P, IDATA, RDATA)
 

Function/Subroutine Documentation

◆ cubrul()

subroutine cubrul ( real*8  F,
real*8, dimension(2,3)  VEC,
real*8, dimension(6)  P,
integer, dimension(1)  IDATA,
real*8, dimension(1)  RDATA 
)
228 C
229 C BASIC CUBATURE RULE PAIR OVER A TRIANGLE
230 C
231 C PARAMETERS
232 C F - EXTERNAL FUNCTION - SEE COMMENTS TO CUBTRI
233 C VEC- MATRIX OF BASE VECTORS AND ORIGIN (INPUT)
234 C P - TRIANGLE DESCRIPTION VECTOR OF DIMENSION 6
235 C P(1) - TRANSFORMED X COORDINATE OF ORIGIN VERTEX(INPUT)
236 C P(2) - TRANSFORMED Y COORDINATE OF ORIGIN VERTEX(INPUT)
237 C P(3) - DISTANCE OF OTHER VERTICES IN THE DIRECTIONS
238 C OF THE BASE VECTORS (INPUT)
239 C P(4) - LESS ACCURATE ESTIMATED INTEGRAL (OUTPUT)
240 C P(5) - MORE ACCURATE ESTIMATED INTEGRAL (OUTPUT)
241 C P(6) - ABS(P(5)-P(4)) (OUTPUT)
242 C
243 C CUBRUL EVALUATES A LINEAR COMBINATION OF BASIC INTEGRATION
244 C RULES HAVING D3 SYMMETRY. THE AREAL*8 COORDINATES PERTAINING TO
245 C THE J-TH RULE ARE STORED IN W(I,J),I=1,2,3. THE CORRESPONDING
246 C WEIGHTS ARE W(4,J) AND W(5,J), WITH W(5,J) BELONGING TO THE
247 C MORE ACCURATE FORMULA. IF W(1,J).EQ.W(2,J), THE INTEGRATION
248 C POINT IS THE CENTROID, ELSE IF W(2,J).EQ.W(3,J), THE EVALUATION
249 C POINTS ARE ON THE MEDIANS. IN BOTH CASES ADVANTAGE IS TAKEN OF
250 C SYMMETRY TO AVOID REPEATING FUNCTION EVALUATIONS.
251 C
252 C THE FOLLOWING REAL*8 VARIABLES ARE USED TO AVOID
253 C UNNECESSARY ROUNDING ERRORS IN FLOATING POINT ADDITION.
254 C THEY MAY BE DECLARED SINGLE PRECISION IF REAL*8 IS
255 C NOT AVAILABLE AND FULL ACCURACY IS NOT NEEDED.
256  implicit none
257  REAL*8 a1, a2, s, sn, dzero, done, dthree, dsix,f,
258  & point5,x,y
259  REAL*8 area, origin(2), p(6), rdata(1), tvec(2,3), vec(2,3),w(5,6)
260  INTEGER idata(1),nquad,i,j,k
261 C
262 C W CONTAINS POINTS AND WEIGHTS OF THE INTEGRATION FORMULAE
263 C NQUAD - NUMBER OF BASIC RULES USED
264 C
265 C THIS PARTICULAR RULE IS THE 19 POINT EXTENSION (DEGREE 8) OF
266 C THE FAMILIAR 7 POINT RULE (DEGREE 5).
267 C
268 C SIGMA=SQRT(7)
269 C PHI=SQRT(15)
270 C W(1,1),W(2,1),W(3,1) = 1/3
271 C W(4,1) = 9/40
272 C W(5,1) = 7137/62720 - 45*SIGMA/1568
273 C W(1,2) = 3/7 + 2*PHI/21
274 C W(2,2),W(3,2) = 2/7 - PHI/21
275 C W(4,2) = 31/80 - PHI/400
276 C W(5,2) = - 9301697/4695040 - 13517313*PHI/23475200
277 C + 764885*SIGMA/939008 + 198763*PHI*SIGMA/939008
278 C W(*,3) = W(*,2) WITH PHI REPLACED BY -PHI
279 C W(1,5) = 4/9 + PHI/9 + SIGMA/9 - SIGMA*PHI/45
280 C W(2,5),W(3,5) = 5/18 - PHI/18 - SIGMA/18 + SIGMA*PHI/90
281 C W(4,5) = 0
282 C W(5,5) = 102791225/59157504 + 23876225*PHI/59157504
283 C - 34500875*SIGMA/59157504 - 9914825*PHI*SIGMA/59157504
284 C W(*,4) = W(*,5) WITH PHI REPLACED BY -PHI
285 C W(1,6) = 4/9 + SIGMA/9
286 C W(2,6) = W(2,4)
287 C W(3,6) = W(2,5)
288 C W(4,6) = 0
289 C W(5,6) = 11075/8064 - 125*SIGMA/288
290 C
291  nquad=6
292  w=reshape((/.3333333333333333333333333d0,
293  & .3333333333333333333333333d0,.3333333333333333333333333d0,.225d0,
294  & .3786109120031468330830822d-1,.7974269853530873223980253d0,
295  & .1012865073234563388009874d0,.1012865073234563388009874d0,
296  & .3778175416344814577870518d0,
297  & .1128612762395489164329420d0,.5971587178976982045911758d-1,
298  & .4701420641051150897704412d0,.4701420641051150897704412d0,
299  & .3971824583655185422129482d0,
300  & .2350720567323520126663380d0,.5357953464498992646629509d0,
301  & .2321023267750503676685246d0,.2321023267750503676685246d0,
302  & 0.d0,.3488144389708976891842461d0,
303  & .9410382782311208665596304d0,
304  & .2948086088443956672018481d-1,.2948086088443956672018481d-1,
305  & 0.d0,.4033280212549620569433320d-1,.7384168123405100656112906d0,
306  & .2321023267750503676685246d0,.2948086088443956672018481d-1,
307  & 0.d0,.2250583347313904927138324d0/),(/5,6/))
308 cparallel
309 cparallel DATA NQUAD /6/, W(1,1), W(2,1), W(3,1)/3*.33333333333333333333333
310 cparallel * 33E0/, W(4,1), W(5,1) /.225E0,.3786109120031468330830822E-1/,
311 cparallel * W(1,2), W(2,2), W(3,2) /.7974269853530873223980253E0,2*
312 cparallel * .1012865073234563388009874E0/, W(4,2), W(5,2)
313 cparallel * /.3778175416344814577870518E0,.1128612762395489164329420E0/,
314 cparallel * W(1,3), W(2,3), W(3,3) /.5971587178976982045911758E-1,2*
315 cparallel * .4701420641051150897704412E0/, W(4,3), W(5,3)
316 cparallel * /.3971824583655185422129482E0,.2350720567323520126663380E0/
317 cparallel DATA W(1,4), W(2,4), W(3,4) /.5357953464498992646629509E0,2*
318 cparallel * .2321023267750503676685246E0/, W(4,4), W(5,4)
319 cparallel * /0.E0,.3488144389708976891842461E0/, W(1,5), W(2,5), W(3,5)
320 cparallel * /.9410382782311208665596304E0,2*.2948086088443956672018481E-1/,
321 cparallel * W(4,5), W(5,5) /0.E0,.4033280212549620569433320E-1/, W(1,6),
322 cparallel * W(2,6), W(3,6) /.7384168123405100656112906E0,
323 cparallel * .2321023267750503676685246E0,.2948086088443956672018481E-1/,
324 cparallel * W(4,6), W(5,6) /0.E0,.2250583347313904927138324E0/
325 cparallel
326 C
327  dzero=0.d0
328  done=1.d0
329  dthree=3.d0
330  dsix=6.d0
331  point5=.5d0
332 cparallel DATA DZERO /0.D0/, DONE /1.D0/, DTHREE /3.D0/, DSIX /6.D0/,
333 cparallel * POINT5 /.5E0/
334 C
335 C SCALE BASE VECTORS AND OBTAIN AREA
336  DO 20 i=1,2
337  origin(i) = vec(i,3) + p(1)*vec(i,1) + p(2)*vec(i,2)
338  DO 10 j=1,2
339  tvec(i,j) = p(3)*vec(i,j)
340  10 CONTINUE
341  20 CONTINUE
342  area = point5*abs(tvec(1,1)*tvec(2,2)-tvec(1,2)*tvec(2,1))
343  a1 = dzero
344  a2 = dzero
345 C
346 C COMPUTE ESTIMATES FOR INTEGRAL AND ERROR
347  DO 40 k=1,nquad
348  x = origin(1) + w(1,k)*tvec(1,1) + w(2,k)*tvec(1,2)
349  y = origin(2) + w(1,k)*tvec(2,1) + w(2,k)*tvec(2,2)
350  s = (f(x,y,idata,rdata))
351  sn = done
352  IF (w(1,k).EQ.w(2,k)) GO TO 30
353  x = origin(1) + w(2,k)*tvec(1,1) + w(1,k)*tvec(1,2)
354  y = origin(2) + w(2,k)*tvec(2,1) + w(1,k)*tvec(2,2)
355  s = s + (f(x,y,idata,rdata))
356  x = origin(1) + w(2,k)*tvec(1,1) + w(3,k)*tvec(1,2)
357  y = origin(2) + w(2,k)*tvec(2,1) + w(3,k)*tvec(2,2)
358  s = s + (f(x,y,idata,rdata))
359  sn = dthree
360  IF (w(2,k).EQ.w(3,k)) GO TO 30
361  x = origin(1) + w(1,k)*tvec(1,1) + w(3,k)*tvec(1,2)
362  y = origin(2) + w(1,k)*tvec(2,1) + w(3,k)*tvec(2,2)
363  s = s + (f(x,y,idata,rdata))
364  x = origin(1) + w(3,k)*tvec(1,1) + w(1,k)*tvec(1,2)
365  y = origin(2) + w(3,k)*tvec(2,1) + w(1,k)*tvec(2,2)
366  s = s + (f(x,y,idata,rdata))
367  x = origin(1) + w(3,k)*tvec(1,1) + w(2,k)*tvec(1,2)
368  y = origin(2) + w(3,k)*tvec(2,1) + w(2,k)*tvec(2,2)
369  s = s + (f(x,y,idata,rdata))
370  sn = dsix
371  30 s = s/sn
372  a1 = a1 + w(4,k)*s
373  a2 = a2 + w(5,k)*s
374  40 CONTINUE
375  p(4) = (a1)*area
376  p(5) = (a2)*area
377  p(6) = abs(p(5)-p(4))
378  RETURN

◆ cubtri()

subroutine cubtri ( external  F,
real*8, dimension(2,3)  T,
real*8  EPS,
integer  MCALLS,
real*8  ANS,
real*8  ERR,
integer  NCALLS,
real*8, dimension(6,nw)  W,
integer  NW,
integer, dimension(1)  IDATA,
real*8, dimension(1)  RDATA,
integer  IER 
)
7 c
8 c changed on June 6th 2011 in order to allow parallelization:
9 c the common statements were commented with cparallel, the data
10 C statements were converted.
11 C
12 C ADAPTIVE CUBATURE OVER A TRIANGLE
13 C
14 C PARAMETERS
15 C F - USER SUPPLIED EXTERNAL FUNCTION OF THE FORM
16 C F(X,Y,IDATA,RDATA)
17 C WHERE X AND Y ARE THE CARTESIAN COORDINATES OF A
18 C POINT IN THE PLANE, AND IDATA AND RDATA ARE INTEGER
19 C AND REAL*8 VECTORS IN WHICH DATA MAY BE PASSED.
20 C T - ARRAY OF DIMENSION (2,3) WHERE T(1,J) AND T(2,J)
21 C ARE THE X AND Y COORDINATES OF THE J-TH VERTEX OF
22 C THE GIVEN TRIANGLE (INPUT)
23 C EPS - REQUIRED TOLERANCE (INPUT). IF THE COMPUTED
24 C INTEGRAL IS BETWEEN-1 AND 1, AN ABSOLUTE ERROR
25 C TEST IS USED, ELSE A RELATIVE ERROR TEST IS USED.
26 C MCALLS- MAXIMUM PERMITTED NUMBER OF CALLS TO F (INPUT)
27 C ANS - ESTIMATE FOR THE VALUE OF THE INTEGRAL OF F OVER
28 C THE GIVEN TRIANGLE (OUTPUT)
29 C ERR - ESTIMATED ABSOLUTE ERROR IN ANS (OUTPUT)
30 C NCALLS- ACTUAL NUMBER OF CALLS TO F (OUTPUT). THIS
31 C PARAMETER MUST BE INITIALIZED TO 0 ON THE FIRST
32 C CALL TO CUBTRI FOR A GIVEN INTEGRAL (INPUT)
33 C W - WORK SPACE. MAY NOT BE DESTROYED BETWEEN CALLS TO
34 C CUBTRI IF RESTARTING IS INTENDED
35 C NW - LENGTH OF WORK SPACE (INPUT).
36 C IF NW .GE. 3*(19+3*MCALLS)/38, TERMINATION DUE TO
37 C FULL WORK SPACE WILL NOT OCCUR.
38 C IER - TERMINATION INDICATOR (OUTPUT)
39 C IER=0 NORMAL TERMINATION, TOLERANCE SATISFIED
40 C IER=1 MAXIMUM NUMBER OF CALLS REACHED
41 C IER=2 WORK SPACE FULL
42 C IER=3 FURTHER SUBDIVISION OF TRIANGLES IMPOSSIBLE
43 C IER=4 NO FURTHER IMPROVEMENT IN ACCURACY IS
44 C POSSIBLE DUE TO ROUNDING ERRORS IN FUNCTION
45 C VALUES
46 C IER=5 NO FURTHER IMPROVEMENT IN ACCURACY IS
47 C POSSIBLE BECAUSE SUBDIVISION DOES NOT
48 C CHANGE THE ESTIMATED INTEGRAL. MACHINE
49 C ACCURACY HAS PROBABLY BEEN REACHED BUT
50 C THE ERROR ESTIMATE IS NOT SHARP ENOUGH.
51 C
52 C CUBTRI IS DESIGNED TO BE CALLED REPEATEDLY WITHOUT WASTING
53 C EARLIER WORK. THE PARAMETER NCALLS IS USED TO INDICATE TO
54 C CUBTRI AT WHAT POINT TO RESTART, AND MUST BE RE-INITIALIZED
55 C TO 0 WHEN A NEW INTEGRAL IS TO BE COMPUTED. AT LEAST ONE OF
56 C THE PARAMETERS EPS, MCALLS AND NW MUST BE CHANGED BETWEEN
57 C CALLS TO CUBTRI, ACCORDING TO THE RETURNED VALUE OF IER. NONE
58 C OF THE OTHER PARAMETERS MAY BE CHANGED IF RESTARTING IS DONE.
59 C IF IER=3 IS ENCOUNTERED, THERE PROBABLY IS A SINGULARITY
60 C SOMEWHERE IN THE REGION. THE ERROR MESSAGE INDICATES THAT
61 C FURTHER SUBDIVISION IS IMPOSSIBLE BECAUSE THE VERTICES OF THE
62 C SMALLER TRIANGLES PRODUCED WILL BEGIN TO COALESCE TO THE
63 C PRECISION OF THE COMPUTER. THIS SITUATION CAN USUALLY BE
64 C RELIEVED BY SPECIFYING THE REGION IN SUCH A WAY THAT THE
65 C SINGULARITY IS LOCATED AT THE THIRD VERTEX OF THE TRIANGLE.
66 C IF IER=4 IS ENCOUNTERED, THE VALUE OF THE INTEGRAL CANNOT BE
67 C IMPROVED ANY FURTHER. THE ONLY EXCEPTION TO THIS OCCURS WHEN A
68 C FUNCTION WITH HIGHLY IRREGULAR BEHAVIOUR IS INTEGRATED (E.G.
69 C FUNCTIONS WITH JUMP DISCONTINUITIES OR VERY HIGHLY OSCILLATORY
70 C FUNCTIONS). IN SUCH A CASE THE USER CAN DISABLE THE ROUNDING
71 C ERROR TEST BY REMOVING THE IF STATEMENT SHORTLY AFTER STATEMENT
72 C NUMBER 70.
73 C
74  implicit none
75  EXTERNAL f,rnderr
76  INTEGER idata(1), ier, mcalls, ncalls, nw,jkp,i,j,k,l,maxc,maxk,
77  & mw,nfe
78  REAL*8 alfa, ans, anskp, area, eps, err, errmax, h, q1, q2, r1,r2,
79  * rdata(1), d(2,4), s(4), t(2,3), vec(2,3), w(6,nw), x(2),zero,
80  & point5,one,rnderr
81 C ACTUAL DIMENSION OF W IS (6,NW/6)
82 C
83  REAL*8 tans, terr, dzero
84 cparallel COMMON /CUBSTA/ TANS, TERR
85 C THIS COMMON IS REQUIRED TO PRESERVE TANS AND TERR BETWEEN CALLS
86 C AND TO SAVE VARIABLES IN FUNCTION RNDERR
87  nfe=19
88  s=(/1.d0,1.d0,1.d0,-1.d0/)
89  d=reshape((/0.0d0,0.0d0,0.0d0,1.0d0,1.0d0,0.0d0,1.0d0,1.0d0/),
90  & (/2,4/))
91  zero=0.d0
92  one=1.d0
93  dzero=0.d0
94  point5=.5d0
95 cparallel DATA NFE /19/, S(1), S(2), S(3), S(4) /3*1E0,-1E0/, D(1,1),
96 cparallel * D(2,1) /0.0,0.0/, D(1,2), D(2,2) /0.0,1.0/, D(1,3), D(2,3)
97 cparallel * /1.0,0.0/, D(1,4), D(2,4) /1.0,1.0/
98 C NFE IS THE NUMBER OF FUNCTION EVALUATIONS PER CALL TO CUBRUL.
99 cparallel DATA ZERO /0.E0/, ONE /1.E0/, DZERO /0.D0/, POINT5 /.5E0/
100 C
101 C CALCULATE DIRECTION VECTORS, AREA AND MAXIMUM NUMBER
102 C OF SUBDIVISIONS THAT MAY BE PERFORMED
103  DO 20 i=1,2
104  vec(i,3) = t(i,3)
105  DO 10 j=1,2
106  vec(i,j) = t(i,j) - t(i,3)
107  10 CONTINUE
108  20 CONTINUE
109  maxc = (mcalls/nfe+3)/4
110  ier = 1
111  maxk = min0(maxc,(nw/6+2)/3)
112  IF (maxc.GT.maxk) ier = 2
113  area = abs(vec(1,1)*vec(2,2)-vec(1,2)*vec(2,1))*point5
114  k = (ncalls/nfe+3)/4
115  mw = 3*(k-1) + 1
116  IF (ncalls.GT.0) GO TO 30
117 C
118 C TEST FOR TRIVIAL CASES
119  tans = dzero
120  terr = dzero
121  IF (area.EQ.zero) GO TO 90
122  IF (mcalls.LT.nfe) GO TO 100
123  IF (nw.LT.6) GO TO 110
124 C
125 C INITIALIZE DATA LIST
126  k = 1
127  mw = 1
128  w(1,1) = zero
129  w(2,1) = zero
130  w(3,1) = one
131  CALL cubrul(f, vec, w(1,1), idata, rdata)
132  tans = w(5,1)
133  terr = w(6,1)
134  ncalls = nfe
135 C
136 C TEST TERMINATION CRITERIA
137  30 ans = tans
138  err = terr
139  IF (err.LT.dmax1(one,abs(ans))*eps) GO TO 90
140  IF (k.EQ.maxk) GO TO 120
141 C
142 C FIND TRIANGLE WITH LARGEST ERROR
143  errmax = zero
144  DO 40 i=1,mw
145  IF (w(6,i).LE.errmax) GO TO 40
146  errmax = w(6,i)
147  j = i
148  40 CONTINUE
149 C
150 C SUBDIVIDE TRIANGLE INTO FOUR SUBTRIANGLES AND UPDATE DATA LIST
151  DO 50 i=1,2
152  x(i) = w(i,j)
153  50 CONTINUE
154  h = w(3,j)*point5
155  IF (rnderr(x(1),h,x(1),h).NE.zero) GO TO 130
156  IF (rnderr(x(2),h,x(2),h).NE.zero) GO TO 130
157  anskp = (tans)
158  tans = tans - (w(5,j))
159  terr = terr - (w(6,j))
160  r1 = w(4,j)
161  r2 = w(5,j)
162  jkp = j
163  q1 = zero
164  q2 = zero
165  DO 70 i=1,4
166  DO 60 l=1,2
167  w(l,j) = x(l) + h*d(l,i)
168  60 CONTINUE
169  w(3,j) = h*s(i)
170  CALL cubrul(f, vec, w(1,j), idata, rdata)
171  q2 = q2 + w(5,j)
172  q1 = q1 + w(4,j)
173  j = mw + i
174  70 CONTINUE
175  alfa = 1e15
176  IF (q2.NE.r2) alfa = abs((q1-r1)/(q2-r2)-one)
177  j = jkp
178  DO 80 i=1,4
179  w(6,j) = w(6,j)/alfa
180  tans = tans + w(5,j)
181  terr = terr + w(6,j)
182  j = mw + i
183  80 CONTINUE
184  mw = mw + 3
185  ncalls = ncalls + 4*nfe
186  k = k + 1
187 C
188 C IF ANSWER IS UNCHANGED, IT CANNOT BE IMPROVED
189  IF (anskp.EQ.(tans)) GO TO 150
190 C
191 C REMOVE THIS IF STATEMENT TO DISABLE ROUNDING ERROR TEST
192  IF (k.GT.3 .AND. abs(q2-r2).GT.abs(q1-r1)) GO TO 140
193  GO TO 30
194 C
195 C EXITS FROM SUBROUTINE
196  90 ier = 0
197  GO TO 120
198  100 ier = 1
199  GO TO 120
200  110 ier = 2
201  120 ans = tans
202  err = terr
203  RETURN
204  130 ier = 3
205  GO TO 120
206  140 ier = 4
207  GO TO 120
208  150 ier = 5
209  GO TO 120
cubrul
subroutine cubrul(F, VEC, P, IDATA, RDATA)
Definition: cubtri.f:228
r1
static double * r1
Definition: filtermain.c:48
rnderr
real *8 function rnderr(X, A, Y, B)
Definition: cubtri.f:212

◆ rnderr()

real*8 function rnderr ( real*8  X,
real*8  A,
real*8  Y,
real*8  B 
)
212 C THIS FUNCTION COMPUTES THE ROUNDING ERROR COMMITTED WHEN THE
213 C SUM X+A IS FORMED. IN THE CALLING PROGRAM, Y MUST BE THE SAME
214 C AS X AND B MUST BE THE SAME AS A. THEY ARE DECLARED AS
215 C DISTINCT VARIABLES IN THIS FUNCTION, AND THE INTERMEDIATE
216 C VARIABLES S AND T ARE PUT INTO COMMON, IN ORDER TO DEFEND
217 C AGAINST THE WELL-MEANING ACTIONS OF SOME OFFICIOUS OPTIMIZING
218 C FORTRAN COMPILERS.
219  implicit none
220  real*8 x,a,y,b,s,t
221 cparallel COMMON /CUBATB/ S, T
222  s = x + a
223  t = s - y
224  rnderr = t - b
225  RETURN
rnderr
real *8 function rnderr(X, A, Y, B)
Definition: cubtri.f:212
Hosted by OpenAircraft.com, (Michigan UAV, LLC)