From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Utilib/src/ALGPT.f | 408 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 408 insertions(+) create mode 100644 Utilib/src/ALGPT.f (limited to 'Utilib/src/ALGPT.f') diff --git a/Utilib/src/ALGPT.f b/Utilib/src/ALGPT.f new file mode 100644 index 0000000..ae46666 --- /dev/null +++ b/Utilib/src/ALGPT.f @@ -0,0 +1,408 @@ +*DECK ALGPT + SUBROUTINE ALGPT(NGPT,XINF,XSUP,ZGKSI,WGKSI) +* +*----------------------------------------------------------------------- +* +*Purpose: +* returns gauss integration points and weights for integration limits +* specified and to the order specified +* +*Copyright: +* Copyright (C) 1991 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): R. Roy and O. Vagner +* +*Reference: +* integral of F(X)DX , from XINF to XSUP is given by summ from i=1 to +* NGPT of F(ZGKSI(I))*WGKSI(I) +* the points ZGKSI and weights WGKSI are generated according +* to the technique described in Handbook of mathematical functions +* M. Abramowitz and I. Stegun, Dover Publication Inc. (1972); +* all points and weights are taken from A.H. Stroud and D. Secrest, +* gaussian quadrature formulas, Prentice-Hall (1966). +* +*Parameters: input +* NGPT number of gauss points +* XINF lower integration limit +* XSUP upper integration limit +* +*Parameters: output +* ZGKSI integration points +* WGKSI integration weights +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NGPT + REAL XINF,XSUP,ZGKSI(NGPT),WGKSI(NGPT) +*---- +* LOCAL VARIABLES +*---- + PARAMETER ( I2= 1, I3= 2, I4= 3, I5= 5, I6= 7, I7=10, I8=13, + > I9=17,I10=21,I11=26,I12=31,I13=37,I14=43,I15=50, + > I16=57,I17=65,I18=73,I19=82,I20=91,I24=101,I28=113, + > I32=127,I64=143,IFN=174) + PARAMETER ( HALF=0.5E0, TWO=2.0E0 ) + REAL RZKSI(IFN),RWKSI(IFN),FW,FC + CHARACTER CERROR*4 + SAVE RZKSI,RWKSI +C N=2 + DATA (RZKSI(I),I=I2,I3-1) / .577350269189626E0/ + DATA (RWKSI(I),I=I2,I3-1) / 1.000000000000000E0/ +C N=3 + DATA (RZKSI(I),I=I3,I4-1) / .774596669241483E0/ + DATA (RWKSI(I),I=I3,I4-1) / .555555555555556E0/ +C N=4 + DATA (RZKSI(I),I=I4,I5-1) / .339981043584856E0, + > .861136311594053E0/ + DATA (RWKSI(I),I=I4,I5-1) / .652145154862546E0, + > .347854845137454E0/ +C N=5 + DATA (RZKSI(I),I=I5,I6-1) / .538469310105683E0, + > .906179845938664E0/ + DATA (RWKSI(I),I=I5,I6-1) / .478628670499366E0, + > .236926885056189E0/ +C N=6 + DATA (RZKSI(I),I=I6,I7-1) / .238619186083197E0, + > .661209386466265E0, + > .932469514203152E0/ + DATA (RWKSI(I),I=I6,I7-1) / .467913934572691E0, + > .360761573048139E0, + > .171324492379170E0/ +C N=7 + DATA (RZKSI(I),I=I7,I8-1) / .405845151377397E0, + > .741531185599394E0, + > .949107912342759E0/ + DATA (RWKSI(I),I=I7,I8-1) / .381830050505119E0, + > .279705391489277E0, + > .129484966168870E0/ +C N=8 + DATA (RZKSI(I),I=I8,I9-1) / .183434642495650E0, + > .525532409916329E0, + > .796666477413627E0, + > .960289856497536E0/ + DATA (RWKSI(I),I=I8,I9-1) / .362683783378362E0, + > .313706645877887E0, + > .222381034453374E0, + > .101228536290376E0/ +C N=9 + DATA (RZKSI(I),I=I9,I10-1) / .324253423403809E0, + > .613371432700590E0, + > .836031107326636E0, + > .968160239507626E0/ + DATA (RWKSI(I),I=I9,I10-1) / .312347077040003E0, + > .260610696402935E0, + > .180648160694857E0, + > .081274388361574E0/ +C N=10 + DATA (RZKSI(I),I=I10,I11-1) / .148874338981631E0, + > .433395394129247E0, + > .679409568299024E0, + > .865063366688985E0, + > .973906528517172E0/ + DATA (RWKSI(I),I=I10,I11-1) / .295524224714753E0, + > .269266719309996E0, + > .219086362515982E0, + > .149451349150581E0, + > .066671344308688E0/ +C N=11 + DATA (RZKSI(I),I=I11,I12-1) / + > .978228658146057E0,.887062599768095E0, + > .730152005574049E0,.519096129206812E0, + > .269543155952345E0/ + DATA (RWKSI(I),I=I11,I12-1) / + > .055668567116174E0,.125580369464905E0, + > .186290210927734E0,.233193764591990E0, + > .262804544510247E0/ +C N=12 + DATA (RZKSI(I),I=I12,I13-1) / + > .125233408511469E0,.367831498998180E0, + > .587317954286617E0,.769902674194305E0, + > .904117256370475E0,.981560634246719E0/ + DATA (RWKSI(I),I=I12,I13-1) / + > .249147045813403E0,.233492536538355E0, + > .203167426723066E0,.160078328543346E0, + > .106939325995318E0,.047175336386512E0/ +C N=13 + DATA (RZKSI(I),I=I13,I14-1) / + > .984183054718588E0,.917598399222978E0, + > .801578090733310E0,.642349339440340E0, + > .448492751036447E0,.230458315955135E0/ + DATA (RWKSI(I),I=I13,I14-1) / + > .040484004765316E0,.092121499837728E0, + > .138873510219787E0,.178145980761946E0, + > .207816047536888E0,.226283180262897E0/ +C N=14 + DATA (RZKSI(I),I=I14,I15-1) / + > .986283808696812E0,.928434883663573E0, + > .827201315069765E0,.687292904811685E0, + > .515248636358154E0,.319112368927890E0, + > .108054948707344E0/ + DATA (RWKSI(I),I=I14,I15-1) / + > .035119460331752E0,.080158087159760E0, + > .121518570687903E0,.157203167158193E0, + > .185538397477938E0,.205198463721296E0, + > .215263853463158E0/ +C N=15 + DATA (RZKSI(I),I=I15,I16-1) / + > .987992518020485E0,.937273392400706E0, + > .848206583410427E0,.724417731360170E0, + > .570972172608539E0,.394151347077563E0, + > .201194093997434E0/ + DATA (RWKSI(I),I=I15,I16-1) / + > .030753241996117E0,.070366047488108E0, + > .107159220467172E0,.139570677926154E0, + > .166269205816994E0,.186161000015562E0, + > .198431485327111E0/ +C N=16 + DATA (RZKSI(I),I=I16,I17-1) / + > .095012509837637E0,.281603550779259E0, + > .458016777657227E0,.617876244402644E0, + > .755404408355003E0,.865631202387832E0, + > .944575023073233E0,.989400934991650E0/ + DATA (RWKSI(I),I=I16,I17-1) / + > .189450610455068E0,.182603415044924E0, + > .169156519395003E0,.149595988816577E0, + > .124628971255534E0,.095158511682493E0, + > .062253523938648E0,.027152459411754E0/ +C N=17 + DATA (RZKSI(I),I=I17,I18-1) / + > .990575475314417E0,.950675521768768E0, + > .880239153726986E0,.781514003896801E0, + > .657671159216691E0,.512690537086477E0, + > .351231763453876E0,.178484181495848E0/ + DATA (RWKSI(I),I=I17,I18-1) / + > .024148302868548E0,.055459529373987E0, + > .085036148317179E0,.111883847193404E0, + > .135136368468525E0,.154045761076810E0, + > .168004102156450E0,.176562705366993E0/ +C N=18 + DATA (RZKSI(I),I=I18,I19-1) / + > .991565168420930E0,.955823949571397E0, + > .892602466497556E0,.803704958972523E0, + > .691687043060353E0,.559770831073947E0, + > .411751161462843E0,.251886225691505E0, + > .084775013041735E0/ + DATA (RWKSI(I),I=I18,I19-1) / + > .021616013526483E0,.049714548894970E0, + > .076425730254889E0,.100942044106287E0, + > .122555206711478E0,.140642914670651E0, + > .154684675126265E0,.164276483745833E0, + > .169142382963143E0/ +C N=19 + DATA (RZKSI(I),I=I19,I20-1) / + > .992406843843584E0,.960208152134830E0, + > .903155903614818E0,.822714656537143E0, + > .720966177335229E0,.600545304661681E0, + > .464570741375961E0,.316564099963630E0, + > .160358645640225E0/ + DATA (RWKSI(I),I=I19,I20-1) / + > .019461788229726E0,.044814226765700E0, + > .069044542737641E0,.091490021622450E0, + > .111566645547334E0,.128753962539336E0, + > .142606702173607E0,.152766042065860E0, + > .158968843393954E0/ +C N=20 + DATA (RZKSI(I),I=I20,I24-1) / + > .076526521133497E0,.227785851141645E0, + > .373706088715420E0,.510867001950827E0, + > .636053680726515E0,.746331906460151E0, + > .839116971822219E0,.912234428251326E0, + > .963971927277914E0,.993128599185095E0/ + DATA (RWKSI(I),I=I20,I24-1) / + > .152753387130726E0,.149172986472604E0, + > .142096109318382E0,.131688638449177E0, + > .118194531961518E0,.101930119817240E0, + > .083276741576705E0,.062672048334109E0, + > .040601429800387E0,.017614007139152E0/ +C N=24 + DATA (RZKSI(I),I=I24,I28-1) / + > .995187219997021E0,.974728555971309E0, + > .938274552002733E0,.886415527004401E0, + > .820001985973903E0,.740124191578554E0, + > .648093651936975E0,.545421471388839E0, + > .433793507626045E0,.315042679696163E0, + > .191118867473616E0,.064056892862605E0/ + DATA (RWKSI(I),I=I24,I28-1) / + > .012341229799987E0,.028531388628934E0, + > .044277438817420E0,.059298584915437E0, + > .073346481411080E0,.086190161531953E0, + > .097618652104114E0,.107444270115966E0, + > .115505668053726E0,.121670472927803E0, + > .125837456346828E0,.127938195346752E0/ +C N=28 + DATA (RZKSI(I),I=I28,I32-1) / + > .996442497573954E0,.981303165370873E0, + > .954259280628938E0,.915633026392132E0, + > .865892522574395E0,.805641370917179E0, + > .735610878013632E0,.656651094038865E0, + > .569720471811402E0,.475874224955118E0, + > .376251516089079E0,.272061627635178E0, + > .164569282133380E0,.055079289884034E0/ + DATA (RWKSI(I),I=I28,I32-1) / + > .009124282593094E0,.021132112592771E0, + > .032901427782304E0,.044272934759004E0, + > .055107345675717E0,.065272923966999E0, + > .074646214234569E0,.083113417228901E0, + > .090571744393033E0,.096930657997930E0, + > .102112967578061E0,.106055765922846E0, + > .108711192258294E0,.110047013016475E0/ +C N=32 + DATA (RZKSI(I),I=I32,I64-1) / + > .048307665687738E0,.144471961582796E0, + > .239287362252137E0,.331868602282128E0, + > .421351276130635E0,.506899908932229E0, + > .587715757240762E0,.663044266930215E0, + > .732182118740290E0,.794483795967942E0, + > .849367613732570E0,.896321155766052E0, + > .934906075937740E0,.964762255587506E0, + > .985611511545268E0,.997263861849482E0/ + DATA (RWKSI(I),I=I32,I64-1) / + > .096540088514728E0,.095638720079275E0, + > .093844399080805E0,.091173878695764E0, + > .087652093004404E0,.083311924226947E0, + > .078193895787070E0,.072345794108849E0, + > .065822222776362E0,.058684093478536E0, + > .050998059262376E0,.042835898022227E0, + > .034273862913021E0,.025392065309262E0, + > .016274394730906E0,.007018610009470E0/ +C N=64 + DATA (RZKSI(I),I=I64,IFN) / + > .024350292663424E0,.072993121787799E0,.121462819296121E0, + > .169644420423993E0,.217423643740007E0,.264687162208767E0, + > .311322871990211E0,.357220158337668E0,.402270157963992E0, + > .446366017253464E0,.489403145707053E0,.531279464019894E0, + > .571895646202634E0,.611155355172393E0,.648965471254657E0, + > .685236313054233E0,.719881850171611E0,.752819907260532E0, + > .783972358943341E0,.813265315122798E0,.840629296252580E0, + > .865999398154093E0,.889315445995114E0,.910522137078503E0, + > .929569172131940E0,.946411374858403E0,.961008799652054E0, + > .973326827789911E0,.983336253884626E0,.991013371476744E0, + > .996340116771955E0,.999305041735772E0/ + DATA (RWKSI(I),I=I64,IFN) / + > .048690957009140E0,.048575467441503E0,.048344762234803E0, + > .047999388596458E0,.047540165714830E0,.046968182816210E0, + > .046284796581314E0,.045491627927418E0,.044590558163757E0, + > .043583724529323E0,.042473515123654E0,.041262563242624E0, + > .039953741132720E0,.038550153178616E0,.037055128540240E0, + > .035472213256882E0,.033805161837142E0,.032057928354852E0, + > .030234657072402E0,.028339672614259E0,.026377469715055E0, + > .024352702568711E0,.022270173808383E0,.020134823153530E0, + > .017951715775697E0,.015726030476025E0,.013463047896719E0, + > .011168139460131E0,.008846759826364E0,.006504457968978E0, + > .004147033260562E0,.001783280721696E0/ +* + IDEP=0 + IFIN=0 + FW=HALF*(XSUP-XINF) + FC=HALF*(XSUP+XINF) + IF( NGPT.EQ. 1 ) THEN + ZGKSI(1)=FC + WGKSI(1)=TWO*FW + RETURN + ELSE IF( NGPT.EQ. 2 ) THEN + IDEP=I2 + IFIN=I3-1 + ELSE IF( NGPT.EQ. 3 ) THEN + IDEP=I3 + IFIN=I4-1 + ELSE IF( NGPT.EQ. 4 ) THEN + IDEP=I4 + IFIN=I5-1 + ELSE IF( NGPT.EQ. 5 ) THEN + IDEP=I5 + IFIN=I6-1 + ELSE IF( NGPT.EQ. 6 ) THEN + IDEP=I6 + IFIN=I7-1 + ELSE IF( NGPT.EQ. 7 ) THEN + IDEP=I7 + IFIN=I8-1 + ELSE IF( NGPT.EQ. 8 ) THEN + IDEP=I8 + IFIN=I9-1 + ELSE IF( NGPT.EQ. 9 ) THEN + IDEP=I9 + IFIN=I10-1 + ELSE IF( NGPT.EQ.10 ) THEN + IDEP=I10 + IFIN=I11-1 + ELSE IF( NGPT.EQ.11 ) THEN + IDEP=I11 + IFIN=I12-1 + ELSE IF( NGPT.EQ.12 ) THEN + IDEP=I12 + IFIN=I13-1 + ELSE IF( NGPT.EQ.13 ) THEN + IDEP=I13 + IFIN=I14-1 + ELSE IF(NGPT.EQ.14 ) THEN + IDEP=I14 + IFIN=I15-1 + ELSE IF(NGPT.EQ.15 ) THEN + IDEP=I15 + IFIN=I16-1 + ELSE IF(NGPT.EQ.16 ) THEN + IDEP=I16 + IFIN=I17-1 + ELSE IF(NGPT.EQ.17 ) THEN + IDEP=I17 + IFIN=I18-1 + ELSE IF(NGPT.EQ.18 ) THEN + IDEP=I18 + IFIN=I19-1 + ELSE IF(NGPT.EQ.19 ) THEN + IDEP=I19 + IFIN=I20-1 + ELSE IF( NGPT.EQ.20 ) THEN + IDEP=I20 + IFIN=I24-1 + ELSE IF( NGPT.EQ.24 ) THEN + IDEP=I24 + IFIN=I28-1 + ELSE IF( NGPT.EQ.28 ) THEN + IDEP=I28 + IFIN=I32-1 + ELSE IF( NGPT.EQ.32 ) THEN + IDEP=I32 + IFIN=I64-1 + ELSE IF( NGPT.EQ.64 ) THEN + IDEP=I64 + IFIN=IFN + ELSE + WRITE(CERROR,'(I4)') NGPT + CALL XABORT('ALGPT:INVALID NBR.OF GAUSS PTS.='//CERROR//' *** + > 1 TO 20,24,28,32,64 PTS. ARE PERMITTED') + ENDIF +C------ +C INITIALIZE ZGKSI AND WGKSI STARTING FROM BOTTOM WITH NEGATIVE +C AND FINISHING WITH POSITIVE +C------ + IUP=1 + IDN=NGPT + DO 100 I=IFIN,IDEP,-1 + ZGKSI(IUP)=-FW*RZKSI(I)+FC + WGKSI(IUP)=FW*RWKSI(I) + ZGKSI(IDN)=FW*RZKSI(I)+FC + WGKSI(IDN)=FW*RWKSI(I) + IUP=IUP+1 + IDN=IDN-1 + 100 CONTINUE +C------ +C FOR ODD NGPT, CENTRAL POINT IS 0.0 AND WEIGHT IS (2.0-TOTAL WEIGHT) +C------ + IF(IUP.EQ.IDN) THEN + ZGKSI(IUP)=FC + DO 110 I=1,IUP-1 + FW=FW-WGKSI(I) + 110 CONTINUE + WGKSI(IUP)=TWO*FW + ENDIF + RETURN + END -- cgit v1.2.3