Numerical inverse isoparametric mapping in 3D FEM

Conprcters & S~lucrrues Vol. 29, No. 4, pp. 6L l-622, primed in Gmat Britain.

NUMERICAL

olw-7949/8% 53.00 + 0.00 6 t988i’ewmonPrcatptc

1988

INVERSE ISOPARAMETRIC IN 3D FEM

MAPPING

V. MURTI, Y. WANG and S. VALLIAPPAN Department of Civil Engineering Materials, The University of New South Wales, Australia (Received 25 June 1987)

Abstract-An

extension of numerical inverse isoparametric mapping [Mm% and Valliappan, Cornput. Struct. 22,101 l-1021 (198611to three dimensional (3D) spaceis presented. In the 3D finite element method (FEM), the determination of the local coordinate of a point of known Cartesian coordinates is not so straight-forward. A direct iterative technique that may require iterations of order 0(N3) is usually employed instead of solving a system of nonlinear equations. This iteration, however, can be improved to be only of order O(Nz) by defining a line in the 3D Cartesian space and iterating along this line. The

details of the method are described and the FORTRAN subroutines for the method proposed are also given.

1. INTRODUCTION

In isoparametric finite element procedure, the evaluation of the elemental matrices (stiffness, mass or damping) should be done numerically using among others the we11 known Gauss-Legendre quadrature [l]. Fortunately, all the required integrations in isoparametric FE procedure can be done in the local (t) or parent space and hence, an inverse isoparametric mapping, (I&f)-‘, defined as (Z&f)-’ s (x,+

5,)

such that

x,,, = N&,)x,,

inverse isoparametric mapping is used. Thus, the latter technique is preferred because of its efficiency and exactness in a sense that isoparametric mapping is used consistently throughout. Implementation of this technique is relatively simple and the method can be realized as a few FORTRAN subroutines, which are given in Appendix 2.

(1)

is never needed. In the above mapping equation, the summation sign has been omitted for clarity, x,, xi are respectively the Cartesian coordinates of a known point (M) and corner nodes, &,,is the corresponding local coordinate of x,, IV$is the shape functions for a 3D element at node-i. The node numbering convention in this paper is shown in Fig. l(a). In certain applications, however, this inverse map ping is inevitably valuable. For example, in remeshing where the initial nodal quantity vectors at the newly created nodes are required [2], these vectors must be interpolated using the values from the previous iterations using numerical inverse isoparametric mapping. In the FE post processing for nodal contouring[3], it was shown that numerical inverse isoparametric mapping is efficient and can be easily implemented. In 3D FE analyses, the stress contours along a certain defined plane of interest are sometimes required, especially in geotechnical fields. This contouring certainly requires the inverse isoparametric mapping to find all the necessary points on the intersection plane, either by direct iterative technique or by numerical inverse isoparametric mapping, to be described in this paper. When a direct iterative technique is used, the iteration required is of order 0(N3). This iteration, however, can be improved to only of order 0 (N*) when numerical 611

(a)

( b)

Fig. 1. Definition of an iterating line PQ.

V. MURTIet al.

612

where

5

I

(W (5’4

82 =(x,),4

-(x&L

and tensor-like notation has been adopted here, where the free index i is used and the dummy indices k, I, m running from 1 to 3 with some restrictions (k # 1,l # n, k # n). xi, (xJi are respectively the coordinates of nodal points and a convenient point P. The constants Ai are defined as:

Fig. 2. Intersection points of two curved lines.

2.

THE METHOD OF INVERSE ISOPARAMETRIC MAPPING

Ai= (xph- (xmh. A similar technique to that presented in [3] is extended here to 3D space. To elucidate the technique clearly, consider an interior point M of coordinate x, = (x,,,, y,, z,,,) of a 3D (finite) element shown in Fig. l(a). Its corresponding local coordinate, c,,, = (e,,,, q,,,, &) to be determined is shown as point M’ in Fig. l(b). Mathematically, the method can be described step by step as shown below. The equation of a line PQ [Fig. l(a)] in 3D space passing through point M can be stated as PQ ~xi-(Xp)=X((Xp)i-(Xe)i),

(2)

where xp, xq are the known coordinates of points P and Q, and x is a constant. The point P is usually taken as one of the convenient nodal points [node 7 shown in Fig. l(a)], of which both the Cartesian and local coordinates are known. Since, point M is also on the line, and only two points are needed to define a straight line, the equation can be rewritten by replacing xp by x,,,, that is, the coordinate of point M. The coordinate of point Q then can be determined by intersecting the line to one of the side planes of the element. The line PQ is in general mapped as a curve line, P’Q’ in the local coordinate space as shown in Fig. l(b). From the isoparametric mapping, the following equations are obtained: x = Ni(5)Xi

(3a)

u = N,(T)&,

(3b)

Thus, for a given plane parallel to one of the parent coordinate planes, 5 = $ = &,(J’= 0 initially), eqn (4) for the line can be rewritten as intersection of two other planes. The corresponding values of the remaining ordinates, (q:),,, (CT),,, can be determined by solving the two equations below:

(74

However, these lines are in general parabolically curved lines (Fig. 2), hence there exist two intersection points. An iterative technique of order (N) can be used to identify such points. Thus, it is obvious that two iterative procedures are required to arrive at the final solution. The outer iterative scheme on tj requires N number of iterations and while keeping the 5 = lj constant, another iterative (the inner iterative scheme) is then required to find the corresponding (r~:)~,([:)~ such that the global iterating points will fall along the iterating line PQ. Instead of adopting the usual bisection technique in the outer iterative procedure, the following iter-

where the summation sign has been omitted for clarity, the vector u is the displacement vector and ui is the nodal displacement vector. Substituting eqn (3a) into (2), the general equation for the curved line P’Q’ can be expressed as the intersection of two planes defined as Ni(5)Yl,i-P1=”

(44

Ni(t

WI

)Y2. i - 82= O9

(6)

Fig. 3. Initial iterating point in the local space.

613

Numerical inverse isoparametric mapping in 3D FEM

i,lp :’

6’

0

:’

,:’

(i

(ii

I

mm._

(iii)

1

Fig. 4. Scheme for the inner iteration.

ative scheme is employed. The initial guess of the iterating point in the local coordinate, 5j= tdj = 0 initially)is set as follows:

the isoparametric mapping whether the Cartesian coordinate corresponding to each pair of the solution lies along the iterating line PQ in the original Cartesian space, such that,

x0*= Ni(t09 (Vi)09(ii)O)xi where &, & are respectively the local coordinates of points P’ and Q’ (Fig. 3), and the proportional constant, 6,, is defined as

6,+

PQ

where

D, = (Ax ;“)o.5

X$E? PQ.

(llb)

If the new coordinate of (to, (v~)~,(t;i)o) lies on the line PQ, the point @ is the initial guess of x,,,. When the solution for the intersection points sought for has not been obtained in the inner loop of iterations, the procedure is repeated by resetting the new values for (vi+ ,)o, (C, ,). using the bisection technique,

(9b) (Vi+

= (Ax;QP)~,’

(9c)

Ax;~ = ((x,Ji - (x,,,)~)*

(9d)

D,

(114

I )O =

o.5*((?i-

I )O +

(rli)O)

(ri+,)O=0.~*((1,-,)O+(ri)O)r

UW

(12b)

until the solution pair is obtained ((v:)~, (cF)o) for a particular value of co. Depending on the shape of the are the Cartesian lengths of lines PM and PQ. Using quadratic curves, the direction of the locus path of to as the initial guess and keeping the value L$,, ((u~)~,(t;i)o) can take one of the paths shown in Fig. 4, constant, the corresponding (v:)~, (I:)o are detersuch that another initial guess may be necessary when mined in the inner iterations using the initial guess path type 2 illustrated in Fig. 4 is encountered during (rt,)o = tlo, ([o)o = Co from eqn (8). Substituting iteration. In general, since the final local coordinate (rt,),, (Lo)ointo eqn (7), the updated sets, which are obtained in the outer iterative procedure cannot be essentially nearer to the final solution of (v~)~,(ci)o, reached in a single iteration, a new iterating point [see are obtainable by solving the following two quadratic eqn (8)] t;, should be selected in a manner illustrated equations (see Fig. 4): in eqn (8) by replacing one of the end points by x$, i.e. the corresponding Cartesian coordinate of point (lOa) (to, (v:)~, (CTlo) (see Fig. 5). The whole procedure is then repeated for each rj until the final solution (lob) CZ = (tj9 (tlF)jv j> satisfies where the coefficients, Ai, Bi, Ci are functions of The complete definitions for ~jrYI,i~Y*,irBI~82~ coefficients A,, Bi, Ci are given in Appendix 1. It is obvious that a multiple solution is in general possible, since both lines are in general curved lines in 3D space. Hence, from a pair of the quadratic equations above, four permutal pairs of solutions are obtained from eqn (10). However, it is easy to check by using

$” =Nk(lj, (~f)j, (CF>j>X,

(W

x,=x;+&

(13b)

for a tolerance &. A typical illustration of iteration involving x7 (i.e. the outer iterative procedure) is schematically shown in Fig. 5. Since the iterations are conducted using the original Cartesian coordinate as

614

V.

MURTIet al.

Fig. 5. Scheme for the outer iteration.

the reference space, it often happens that the magnitude of the intended coordinate (of point M) is vastly greater than the length of the element or vice versa; hence, different tolerance limits must be employed in the outer and inner iterative procedures. The corrected tolerance limits may be employed by firstly defining, ai = 0.5*(x? - xy},

Wt

where xyBx,xy are respectively the maximum and minimum nodal coordinate components. The two limiting tolerances can then interchangeably be defined as follows:

A point M is an interior point if its projected coordinate on each plane is also an interior of the clipping plane on that plane. An axis tra~sfo~ation may be necessary in this case to obtain a suitable orientation or alignment with the element concerned. Since a general 3D (finite) element may have curved planes for all its six sides, their projected planes on the auxiliary planes will have necessarily curved sides. The dete~ination of the clipping planes then becomes unnecessarily tedious. Thus, an alternative technique can be devised easily by using the basic property of set theory. From elementary set theory, a point x,_~ = (x, y) belongs to a union of sets A and B (i.e. the clipping plane), if it belongs to set A or set B. Mathematically, AUB=A$B@(AnB) x,_,e(A

UB) = (x,_,,tcA)

aN’dN,l

or (x,_~E~?).

(1Sb)

Hence, if Ex_,,, ” indicates the two projected planes (n = 1,2) of a 3D (finite) element E on the auxilliary plane x -y, a point M is an interior point if (x:-~%-,J

and

(~!LEJL,,)

and

(J$‘_~EE~_~,~)= true,

where the condition 6N. 2 =

(154

(16a)

in each bracket is

W4

aN can be of any appropriate magnitude, In the subroutines given, a value of 70 has been adopted for a,, but this value can be arbitrarily modified deRending on the accuracy sought. This completes the whole procedure of the iterative technique of the numerical inverse isoparametric mapping. Since two sets of iterations are required in this procedure, the rn~irn~ number of iterations for a given tolerance, S,, is of order O(P) compared to the direct iterative technique of order O(N3). The procedure of this technique is given as FORTRAN subroutines in Appendix 2.

where

3. LOCATING AN INTERIOR POINT IN A 3D FINITE ELEMENT

A similar technique for locating an interior point to that presented in [3] can also be extended to a 3D case. In a two dimensional element, an interior point can be located quite simply by traversing all elemental sides and t~nsfo~~ng the local axis accordingly as illustrated in Fig. 6. The details of this method already given in [3] will not be repeated here. However, the method is simple and the equation can be constructed in a straightforward manner. In a 3D space, an interior point can also be assessed by inspecting all the three auxiliary projection planes, that is the x-y, x-z, and y-z planes.

or

(x;-~E%,.,z),

(16b)

and XT_? indicates the coordinate of point M on the x - y auxilliary projection plane. In a brick 3D element [Fig. 7(a)], the clipping plane may coincide with the projected plane, Ex_,,+ or may consist of a union of only two simple rectangles, such that the inspection of interiority becomes trivial. However, the above method is robustly valid only for general concave outward elements (brick elements included) [Figs 7(a) and (b)]. Certain exterior points as seen in Fig. 7(c) surrounding a general 3D concave inward element will not be detected using the above method. Hence, a more rigorous check for such elements is required. When concave inward elements as shown in Fig. 7(c) are used in the analysis or possibly generated by a remeshingladaptive routine, the following method may be employed. Considering an exterior point 0 shown in Fig. 8(a), a plane parallel to the zaxis and containing the point 0 can be easily constructed as the 3D clipping plane (plane C). When the key points defining this clipping plane have been identified, the intersection shape of the plane C and the element can be reproduced as typically shown in Fig. 8(b). From this point onwards, the inspection technique presented ~h~aticalIy in Fig. 6 or described in [3] can then be employed to assess the interiority of the

615

Numerical inverse isoparametric mapping in 3D FEM 5

(b)

08

Fig. 6. Determination of interiority of a point M.

point 0 with respect to the element [Fig. 8(b)]. The only tedious step of this method is finding the intersection plane mentioned earlier. However, in spite of the tedium involved in analytical techniques to obtain this plane, hardware oriented algorithms with the aid of computers with graphics capability are available, making this step trivially simple. One possible technique using the projected plane to find one side of the intersection plane is shown Fig. 9.

4. APPLICATIONIN STRESS CONTOURING In geomechanics problems the stress levels on a certain defined plane are often required for analysis, especially in analyses of foundations of large arch dams where the existence of a plane of weak discontinuity cannot be readily assured. In such analysis, the stress levels on a possibly potential plane (Fig. 10)

V. Muon er al.

Fig. 7. Brick, concave outward and inward solid elements.

should be inspected from a continuum analysis by drawing contours of stresses on that plane. Since the plane can be arbitrarily oriented in any direction in a 3D space, and since the Cartesian space is usually used as a reference space, only Cartesian coordinates

Fig. 9. Determination of solid intersection surfaces.

of the key points of the intended plane are available. However, the stress levels at the points on the plane can be determined only by first finding their corre-

(01

(b)

Fig. 8. Determination of interiority of a point in three dimensions.

Numerical inverse isoparametric mapping in 3D FEM 5. DISCUSSIONS

sponding local coordinates in the parent space. The numerical inverse isoparametric mapping is an appropriate technique not only for this step, but also for locating the intended plane of an arbitrary direction in 3D space where stress levels are required. Once the local coordinate is determined, the stress level at that point can be calculated using the following:

Acknowledgements-The work described in this paper forms a part of research programme on wave propagation and fracture mechanics problems being carried out in the Department of Civil Engineering Materials, The University of New South Wales. The authors wish to acknowledge the financial support received through the Australian Research Grants Scheme. REFERENCES

K. J. Bathe and E. Wilson, Numerical Method in Finite Element Analysis. Prentice-Hall, Englewood Cliffs, NJ (1978). S. Valliappan and V. Murti, Automatic remeshing technique in quasi-static and dynamic crack propagation. NUMETA-85, Swansea, Proc. (Edited by J. Middleton and G. N. Pande), pp. 107-116. Balkema, Rotterdam (1985). V. Murti and S. Valliappan, Numerical inverse isoparametric mapping in remeshing and nodal quantity contouring. Comput. Struct. 22, 1011-1021 (1986).

(17)

where a,,,, &,, are respectively the stress vector at a point M and the local coordinate of point M obtained by using the numerical inverse isoparametric mapping, 5: = (rj, (r7t)j([T)j). The matrices [Dl, and [B], are the stress-strain and strain-displacement relationships, respectively, and II; are the nodal displacement vectors for an element e.

APPENDIX

AND CONCLUSIONS

In a 3D FE analysis where only brick elements are used, the technique may unnecessarily become inefficient, since as mentioned earlier in the paper the procedure is considerably simplified and the solution becomes trivially simple. Since generality is a feature of the method presented, it is applicable to any shape of 3D isoparametric elements. In remeshing, nonlinear or large displacement analyses, this technique is inevitably valuable. It is also particularly useful in 3D FE post processing, especially involving curved 3D elements, such as in arch dams with rock foundation of weak discontinuity. The method presented is general, efficient, systematic and exact in a sense that isoparametric mapping has been used throughout FE procedure. Implementation is straightforward and its application to stress contouring in a 3D FE anaysis has been schematically given in [3]. It is also believed that salient applications of this technique to other aspects of FE analysis are possible.

Fig. 10. Arch dam coarse FE discretization.

0, = PleP(Lmc,

617

I: COEFFICIENTS

A, B

AND c

Table Al. Explicit expressions for coefficients A, B, and C Node No. (i)

ai

bi

Comer nodes

(1 + %)(I + 0Yt,i 8

(tt, + r, - 1)(1 + L)(l + rl,)&Yt., 8

2,6, 14, 18

-(J + %)(I + 0Yk.i

0

4

9,10, II,12

0

4,8, 16,20

0

<,=L’t9

~~=~i’tl,

(1 - r1*w + U&Yk,i 4 (1 - r21u + rlJ5Yk.i 4

L=Ci.L A,= i q, 4~ j-l

i /=I

ci (n,+4,,-2)(1

+i,)(J 8

+rl,)Yk,i

(1 + &)(l + C)Yr,; 4

(1 - r1*)(1+ CJYk,i 4 (1 - 12)(1 + rL)Yh,i 4

bp Ck= f Cj-Bk. j-l

The expressions for /JI and yk,i are already given in eqn (5) of the present paper. The above expression is given for k = 1; the expressions related to indices of k equal to 2 and 3 can be obtained in a straightforward manner by cyclical permutation of the r,, ranand c, parameters.

618

V. MURTI et al. APPENDIX 2: FORTRAN SUBROUTINES

c

c------computetheparameters of the twocurved swybce equations c i.e. (rxl,rx2,cOl,cO2) subrouttnea2r3d(.vn,cno,nenxlm) c______________-__________________________________________---------c c ptWp0se: 60 ~a~~rxCrrxl~~,cm~l~.cm(2~~~3J,cn0p(l~,~~~2J,cn0p(3~ * Routineto@d the local coordbtateof an interiorpoint c cno.lZ.3.cOl .cO2,nen) got0 90 &nown Cartesian coordbtate c 1----

c btputparameters: xm(3)--hnownCartesiancoordbtateofan interiorpoint- M c cno(3,nen)-Cartesian coordimtes OftmaWpobtts E nen-number of nodes in an element c c output parameters: tim(3)--local coordinateofpoint M c

c

70 call rxc(rxlr~,cm(2),cm(3)~~lJ,~~2J,c~~3J,~o~l) * .cno23.l.cOl.cO2.nenJ got0 90 c 80 call rxc(rxl,r*2.cm(3).cm(l)cm(2),~o~3).cnop(l),~o~2) * .cno.3.12~-Ol,dlz,nenJ

c c------calculateintersection coordbtateofintermidiate lines, MC & PH c routine calls: c rxc,itertl ,itert2 c 90 do 100 ill.3 c ________________ ______________________________________ _~~_~~~~~~~~~~~ cm(i)=cno(i.7)-cno(i.IS) c 100 continue implicit rear8 (a-h,o-z) c aa=cm(l)+cm(l)+cm(2)+cm(2)+cm(3)*cm(3J dbnensionc?d3,nenJ.vnf3),xim(3),cnop(6).cnoq(6),cm(6), * bb=cm(lJ+(xm(lJ-cnopO)+cm(2)“(xm(2)-cnop(2)J+ rx~f~~~~(~J,gsp(3)~~~20)~~~(6)d~~3), * + cm(J)+(.un(3)-cnop(3)) cnomaxf3).cnomin(3),cll(3),cg1(3) tt=bblaa c c c---initializations do 110 i=l,3 c cm(iJ=cnop(i)+tt%n(iJ nn=o 110 contitute n2=0 c tt3=1 c------estimate local coordinatesofpoints G and M hh=o c c dpg=dsq~(cnop(l)-cm(l))~(cnnp(l)-cm(l))+(c~~2)~~2))~ c------checkifpoint M and point P (stam’ngpoint) coincide * c (cnop(2)-cm(2JJ+(c~~3~~3J~(c~~3~~3JJJ dph=dsqrt((cnop(lJ-cno(l,l5))+(cnop(lJ-cno(l,l5))+ ~#xm(l).eq.cno(l,7J.and.xm(2J.eq.cnoi2,7J.and.and3J.eq. * (cnop(2)-cn0(2,15))*(cnop(2)-cn0(2.15))+ Cno(3.7)) got0 190 * (cnop(3J-cno(3,15JJ*(cnop(3J-cno(3~5JJJ c c----find he maximtonand minimumvaluesofcoordinates of nodalpobtts dgm=dwMmdlJ-cm(l JJ+fxndl)_cm(lJJ+(xm(2J-cm(2JJ* *c of element (a(z)-cm(21J+(N3J-n(3JJ*(N3J&3JN if(dgm.eq.O.0)dgm=l .O c c do 10 1=1,3 do 120 i=l3 cnomax(iJ=cno(i,lJ cml(i)=cnop(i+3)-(cnop(i+3)+cnop(i+3)J*dpghfph cnombt(i)=cno(iJ) gsp(i)=cml(i)+(*m(iJ-cm(i))*2.O/(dgm*(c~iJ10 contbtue * c CnominoJJ 120 continue do 20 i=l ,nen c do20]=13 gsp(ii)=l.O iflcnomax(l).lt.cno(j,iJJcnomaNjJ=cnofj.i) c flcnomin(j).gt.cno(l.i)) cnomin(jJ=cno(j,iJ c----calculate the coordinateofpoint Q 20 continue c c call itertl(n2,J,ii,rxl ,r.&Ol ,cO2,8sp(l),gsp(2),gsp(3). c------calculatethe limitedvaluesfor convergence * Cll) c do 25 i=l.3 call nn20(gsp,shap,nenJ cil(i)=l .Oe-6*(dabs(cnomax(i))+dabsYcnomin(i)))l2.0 c cgl(i)=7OSPcll(i) do 140 j=l3 25 con&me Cnog(/+3)=8sp(l) c cnoq(j)=O.O di=c8l(l)+c8l(lJ+c8I(2)“cgl(2)+c8l(3)+cgl) do 140 i=l,nen c cnoqfj)=cnoq(j)+shap(iJ*cno(j.iJ c------chooseintial plane for iteration 140 continue c c do4Oh13 c---cakuhue lengths ofPh4 and PQ cnop(t)=cno(i,7) c cmfi)=cnop(i)-am(i) 150 dpm=(cnop(l)-.am(l))*(cnop(l)-xm(lJ)+(cnop(2J-xm(2))* cnop(i+3)=-1.0 * fcnop(2)-~2J)+(cnop(3)-~3J)*(cnop(3)-~3)) 40 continue dpq=(cnop(l)-cnoq(l))*(cnop(lJ-cnoq(l))+(cnop(2)-cnoq(2J)+ c * (cnop(2)-cno9(2J)+(cMP(J)-c~3))*(~~~3)-~~(3)) ccm=cm(l) c iid c-----calculate the initial position for iteration c c do 50 i=l3 if(dpq.gt.dpm.and.dpq.le.0.1)goto 160 fldabs(cm(i)J.le.dabs(ccm))goto 50 dd=dsqrt(dpmldPqJ ccm=cm(iJ xsin2=cnop(4)+dP(cnoq(44(4tcnopo) ti=i etM2=cMp(5)+dd+(~(5)~~~5)) 50 contbme csin2=cnop(6)+dd’@(cnoqf6)-cnop(6)) c got0 170 goto@O,7Ot9OJ. ii

.,

Numerical

inverse isoparametric

mapping in 3D FEM

619

c

c

160 xsin2=(cno~4)+cnoq(4))Iz.O etan2=(ctwp(5)+cnog(5)HZ.O csin2=(cnop(6)+cnoq(6))/2.0 C ~---fist

iterative

precedure

return

end subroutine rxc(rxl~~.nnl~,nnt,cnl,cnl~~~,cno, * kl&2,k3,cOl,cO2nen) c______ ______________ - ____~~~~~~~~~~~~~~~~.~~.~.~~~~~~~~~~~~~~~~~~~~

c 170 call ite~l(~~,ii,rxl,rx2,cOl~2wi~,eum2,c~i~, * cll) c c--c

second iterative procedure

call itert2(~rm,Mdpmsl~dcno~~p,gsp,cnop,c~, + ~rsin2,etan2,csin2dcnocm,nnl,cgl,kcy) flkey.eq.0) got0 150 c do 180 Cl3 xim(i)=cm(i+3) 180 continue got0 200 c c---when point M andpoint P coincide c 190 do 195 i=l3 xbn(i)=-1.0 195 continue c 200 return end subroutine nn2O(gsp,shap,nen) ~_____________________________ ___-_______ _________~~~~~~~~~~~~~~ c parpose : c Shape function routbu for ZO-nodedelements &_I_______ - ____-----___ _____________ _____________ ______ ___________ __ c inpw parameters: c gsp(3)-local coordinates c nen-number ofnodes in an element c outputparaouters: c shap(nen)--values of shape functions at nodes

Erouthe called by: c

tx2r3d,itert2

c

implicit rear8 (a-h,o-s) dimension gsp(3).shap(nen)da(3)db(J) c do 10 i=1,3 da(i)=(I .O-gsp(i))l2.0 db(i)=(l .O+gsp(i))l2.0 10 continue

c pwpose: Routine to calculate parameters ofthetwoctoved swfies c c i.e. (rxl,rx2&l,c62). ~________~__~~~___~_________~---------~~~~~~~~~~~~---------------~~~~~ c input parameters: cm1 ,cm2,cm3--direction cosines oflinePM c c en1 ,cn2,cn3--Cartesian coordinates ofpoint P cno(3,nen)--Cartesian coordinates ofnodal points c c kl,k2.k3--coodinate number I-X, 2-Y. 3-2 ~_______~---~.~~~___________~--------~~~~~~~~~~~~~-----------~-~~~~~~~ c ourput parameters: rxl ,rx2,cOl,E02--parameters of two cwved sta$ace equations c c--- _______I___ --- __~~I~_~_~________~~~~~~~~~~~~~~_____~ c routine called by: txZr3d c c______ __-___ ________________-----__ ___ ___________----__ ___ ___________ c hplicit real*8 (a-ho-s) dimension rxl(l),rx2(l),cno(3,1) c i~dabs(cml).lt.l.Oe-5) goto 50 c c------if line PM notpararell to plane (cnl=const.) c ak=cm2icml Co1=ctL?-ak%nl c do 10 i=l ,nen rx1(i)=cno(k2,i)-a&no(kl,i) 10 continue c ak=cm3icml cQ2=ct&ak+cnl c do 20 i=l ,nen rx2(i)=cno(~,i)-akcno(kl,i) 20 contbtue got0 70 c c------if line PM is pararell to plane (cnl =const.) c 50 cOl=cnl cO2=en2 c do 60 i=l .nen rxl(i)=cno(kl,i) rx2(i)=cno(kl,i) 60 continue c

a=l.O-gsp(l)+gsp(l) b=l .O-gsp(2)+gsp(2) c=l .o-gsp(f~gsp(3) shap(2)=a*da(2)*db(3) shap(4)=db(l)%(2)*c shap(6)=a*da(2)%ia(3) shap(8)=da(l)+da(2)“~ shap(9)=da(l Pb*db(3) shap(lO)=db(lPb+db(3) shap(ll)=db(l)‘b*&(3) shap(l2)=da(l)+b”da(3) shap(l4)=a*db(2)*db(3) shap(l6)=db(lPdb(2)% shap(l8)=a+db(2)%ia(3) shap(2O)=da(lPdb(2)*c

70 return end subroutine pmt(ii,rx,&emn,csin,a,b,c) c_______ _____________ ________ __________ __________ ____________ _________ c purpose : I.7 Routine to calculate parameters (ARC) . c_______.--__________________-________________________________________ c input parameters: c ii-operative plane offist iterative procedure c rx,cO-parameters ofa ctaved swface equation c etan,csin--local odinates ~-~~~~__---------~~~~~~____~------~~~~--______~~~~~_______________ c OUQUtpOramnrrs: c ab,c-parameters of the quadratic equation ~~-~~~~_----------~~~~~~___-------~-____~~~~~~__~~~-___~~~~~~~~~__ c routine calls: c nouns ~~~~_~~_------~-~~~~~~__~____~---~~~~-___________~_________~~~_______ c routine called by: c local ~~-~-~~~.---------------~~~~~.---------~~~~~~~~~~~--~~~~~~~~~~~~~~~~~_

620

V. MIJRRef al. c

200 do 250 i=l20 r.xfi+2O)=rx(mapl(i)) 250 continue got0 400

bnplicir rear8 (a-ho-z) dimension rx(1) a=0 b=O c=-CO

c 300 do 350 i=I 20 rx(i+2O)=rx(map2(i)) 350 continue

call ntrans(rxJi) E

400 return end

subroutine root(a,b,c,rlt) ~“--~--_-~~~---~_________~_-_~-____---~~~~~~”~~~~_~~~~______~~~_______ c purpose: c Routine to solve quadratic eqaation ~_---.___-----____________________----_-~~~~.-___________~~~~~_______ c input parameters: c a,b,c-parametes of quadratic equation ~_____.____________________________---_____-*________________-_______ c output parameters: c r&-one valid root of quatbatic equation c routine called by:

c local ~_**~-._______-____~~_____________---_~-~~~**~___________~__~_”~*~~~ c implicit real*8 (a-h,o-z) c C

~dab~a~.It.l.Oe-5~ goto 510 d~=b~~4.O~a~~ ijldisc.lt.6) gore 510 c disc=dsqrt(disc) C

r~tl=~-~dis~~(2.0~~~ rlt2=(-b+discf/(2XPa) c=.c+yyl*(rx(25)+r~27)~+yy2*(rx(2l)+rxo)+yy3 ifldabs(rltl).ledabs(rlt2)) goto 500 * +yy4*(rx(33)+rx(35)) r=rltl ~c+zzl*(rx(24}+rx(28))+zzZ*lrx136)+rx(4011+lz3c+ rltl =rit2 * ~30~)+224~{r~3l)+r~32)~ ruZ=r c=c+xxl*(rx(26)+rxl26))+x3c/rxo+rx(22~+rx(22j}+~*fr~38)+r~38~~ c * +xx4*(rx(34)+rx(34)) 500 rlt=rItl got0 PO0 return c end SIO t~~b~b).~t.i.Oe-5~ goto 520 &=-c/b subroutine ntruns(rx,ii) &--.--------- .___. --_ ~~~~~~~~.~~~~~~~~~~~~~~~__~~.~~.__ got0 900 c c purpose: c 520 write~6,2000) Routine to transform the local coordinate system w.r.t c the operative p&me of iteration (xsinetan and const. csinj 2000 format(lx.’ NO VALID RESULT - PR~~~~~5D ‘I c _--_-*__---- --______ -_~~~~~~~~~_~_.I.__~_~I__ _____ *._*___--______ stop c c input parameters: c 900 return ii--operative plane for first iterative procedure ~~~---~-_-------~~~~~~-~--------_~_~~~~.__~~~~~~______~.___.~~~~___~~~ ~~routme locol(nl,~n2rxI,rk?,car,c02Sf,rI~l,~~2~,dI~, c upakte parameters: + c r~4O~--~a~ters before and after rotation ~coordi~te srhi,srh2&.n3 jj) ~~~~~~__..~~~--__-_~.~~~~~~~--__~-~~-~~~~~~~~__~~~_____~~*~...~~~~_ C system ~“~----~______-.~~~~~~-~__________~~~~~.~~~~~~~~____~~..__._~~~~____~~ c purpose.. c Routine to calcuiate the new iocal coordinates during c routine called by: c jirsr iterative procedure . c pmt ~.~~*____**~~*__-______~-~--__________~~~*~.___-_*___~~~~*_______ c_** -_-________._ _***_-__ _________ ____*_________ _____._.______ _______ c input parameters: c c ml m2--2--local axes number (12) , (23) or (3J) bnplicit real*8 (a-h,o-z) c rxl ,rx2,cOl ,c#2--parameters oftwo curved sut$zcc equations dimension rx(l),mapl(2Ohmap2(20) c n2,n3--number @iterations reached c ~~~---_*~~-----_-----~~~~--~--------__-~~~~.~.~~~~~~_____~~~.~~~~~~ akm mapliS,liJ7,18,l9,12,7~6,4,162OtObj,IOJJ,14,l3,9~~1 c update parameters: data mapZ119~0~39J~,7,12,18,14f,6~7,16,15,10~,4J,Iit sl,rl ,hl,s2,r2&2--local coodinates on two cuwed su&ees c c c dl,a%direction of iterative increment goto(lOO3tXA300), ii ~----.~~~~-----------~~~~~~----------__~~~*.~~.~~~~~~___~~..”~..~~~~~~ C c routine calls: 100 do 150 i=lfO c pmtroot r~2O~i~=r~i~ ~~~~*.___~~~*-~__~~.~*~~~~~~~.~~~--~~~~~~~~~~~~~____~~.~~.~~~ 150 continue c routine calied by: got0 400 c itertl I

__.--____

-

NumericA

*

cii(3~,or~b~e~-c6onr~.sr,clU2~~ got0 200

bnpllcit

real%?

621

inverse isoparametric mapping in 3D FEM

ia-h,o-z)

goto I10

c

dimension rxl(l),rx2(1)

I70 xsinl =xsin2 ctani =etan2 ~S~~=~~ x&t2=xsbl3 etan2=eran3 csbG!=csin3 ~~~~~~-~inll,gt.cll~i~.or.~bslcsi~~sini~.gt. * cl((3~.or.~b~eta~-etani).gt.c11(2)~ goto II0 got0 200 c

180 xsi~=(mrinz4~~~~~i~)l3.0 etan2=(etan2+etanl*etan3)l3ranj)/3.0 csin2=(csin2+csinl+csin3)l3.0 c 2W return end subroutine itert2(nnnen,kk,dpmdpq,dl,cno,shap,gsp,cnop, return * cnoqxmxsin2,etan2~sin2&no,cm,cml, end * cgtJrevt subroutine irrrtl(nt~,it,rxfxx2,c0I,c02~,eta~,csin2. ~-~-“-“~~~~~~---____~~~~.~_-~_____.-_---~~~~___~~~~~__~______.~~~~_ * Cli) c purpose: ~~--____~~~~~~__---~--.~*--___-----------------~~~.~~~~~~~.~~~”~_~~~~ c Routine for second iterative procedure cpurpose: ~---~~*--.----*~-~-**--~~.~--~-~.-----^-~~~~~~~~~__~~~~~~.~~c Ro~‘~~r~~ iterative procedure . c inputparamerers: ~-~__~~~-___~------~~-~~__--___~~-*-----~~~~~~~~~~~*~~~~~”___~*c dpqdpm-distances of PQ and PM c inputparameters: E cqi(3W-vaiues ofiterattve convergence limits c ii-operative plane for first iterative procedure c cn~6),cno~6)~~3~-~~rdi~tes of point PQ and hi c rxl sx2,col ,cO2-+rameters of twocurved swface equations c xsin2,etan2,csin2-local coordbuttes of iterative point c cll(3)-values of convergence limits dcno(3)-direction of iterative increment ~.-.-___.~~~~~__----~----____~----~.----~~-~~~~~*~~~~~.~~~~~_~~~~ c c cml(6J--coordmates qfpoint M in last step of iteration c update parameters: c ~~3,ne~)--Ca~re~an roordinates of ~1~~ E &&-mtmber ~iteratio~ reached c nen-number of nodes in an element c xsin2,etan2,&&-local coordinates of iterativepoint ~-----.-----_----“~~~.*---_~-~~-___~----~~_-~~~~~_~~~~___~~______~~~__ ~...-___.~~_____~---~-~-~_~~__--___-~.--~~___~~~~~~.~~~~~~~.~~_____~__ c output parameters: c rot&e calis: c cm(6~-coordinates of point hi in this step of iteration e &co1 c W-number of variations ofdirection of iterativetncr. ~.--____.~~~.*__~--~-~-~~~~~~I--__~~.--..-~~~~~~._~~~~~~~~~_____~~~. c nn-number ofseconditemtion reached c routine called by: c hey--flagfor convergence key=1 converged C &&3d ~---~~~---.~~_~----““.~~~~._.______~.-~~~_~~_____~~__________..___“~ ~.~._____~~_____~-----_~~*~~-_-__~~*---*~*~~~~~~_~~*~**~~_____*_~~_ c rout&e calied by: C c tx2r3d bnplictt reai%S (a+-r) ~--..~~-----_~~--~.-~~~~.*-__-__--~-*~~~.~_____~~~_______________” dbnension rxi(i~?~~l~~clI~~~ c implicit real% (a-ho-z) n3wl24f dime&on cno(3.f)~~~l),gs~ssp(l).cnop(l),cnoq(l)), J=O * ~i~l.cmliI~~no(l),cgl(i~ c 110 ?&?=?l24i x&3=xsin2 etanl=etd CSitd-CSid e

goto(l2O,I3OJ40)$ C

I20 call locayz$,rxl,rrz,co~,Eo2xrin2sMnz,csinzx,eta~, * csin2detan~~s~n,eta~,cs~~~~~ jj) goto I30 c 130 cat1 t0ecrlrI9rx2,rxi.~2sot.csin?xsin2,~~~~i~,et~, * csinzdtsindcsmxJtn3~st~~~,~~~ got0 I30 c 140 call local(l2fxl~~,cOI,dl2,csi~xs~~,eranZ,csi~~i~~ + eranz~~,dctan~i~,eta~,~~~~) c 150 iflfi.ge.2) got0 lB0 &j.ge.l) got0 160 got0 170 c 160

dwl24I

U=O xsb&?=(xsbt34.&nl~lz.O eta~~e~iem~)~2.0

csin2=(csin3+csin~1/2.o ~~~i~-rFinl).gt.~R(~~.or.~bsics~~~sin~).gr.

nn=nn+s gsp(l)=xsin2 gsp(2)=etan2 gsp(3)=csin2 c----catculate shape functions c call n~O(gsp~~,nen~ E c----caIcuiate Cartesian coordinates of &erativepoint G c do26Oj=l,3 cm(j)=O.O do 250 i=l,nen C~j~~cm(i)+shapriI*cno(ij) 230 eontbtue 260 continue c

622

V. MURTIet al. i~~.gt.l~~~c~h.ne.dcno(l)) dcnof iJ=dcnoh 310 continue

E

c

i&ipq.gt.dpm.and.dpg.gr.dpmJgot0280 iAdW.le.dpm.and~rn.~.~rn) goto 280

c----move

pointQ to poiitt G

c do 270 i=I .3 cnop(iJ=cmfi) CnOp(f+3)=gSp(iJ

270

conthue

gOtO 300 c c-----move point P to point G

c

280 do 290 i--1.3 cnoq(i)=ndiJ cnoq~i+3)=gsp(iJ 290 continue c c------testif convergence is reached c 300 ~~bscnncl)-~(~)).Ir.cgl(lJ.and.dabs(c~2)-~2JJ.~. * C81(2).and.&bsfnd3)-~~3JJ.l,cgU3)) gOtO 350 ijLImg.It.dlJgOtO 350 iflnn.le.1) got0 305 iffdab~c~l)al(I)).I.cgyll.and.dab~~(2)-cm~(2JJ.~t * cg~(2j.a~.&b~cmC3J-~~~3JJ.~t.cglo) got0 350 c c----cakulote direction of iterativeincrement c 305 do 310 is13 dcnoh=cm(iJ-xm(i) dcnoh=dsignf 1 .O,dcnohJ

+&+I

c c------stopiteratingif number of direction variationeu c iflkk.ge.2)8OtO 330 do 320 i-I,3 cml(i)=cmfiJ cml(i+3)=gsp(i) 320 continue c c----futcrl iterationhas not reached c key-o gOtO 400 C

330 do 340 i-1 $ cm~iJ=(cmfi)+cml(iJ.Jl2.O cm(i+3)=(cml(i+3)+gsp(iJ)l2.0 340 continue key=I C SOtO &o C

350 key=1 do360 i=1,3 cm(i+3J=gsp(iJ 360 continue c 400 return end