Technical Manual - Brunel University London
Contents
1. 2 k 1 k 2 i gt aapi l 1 This leads to the following linear system 2 k 1 2 mil i gt aamj D I 1 If we define bj m 15 then we have aa a 05040 For computational efficiency we later extend the coefficients a with respect to y to a basis of Qx 1 x 2 K which eliminates the need to distinguish several special cases in the evaluation of amp 10 2 2 Raviart Thomas elements RT on triangles P K all polynomials of total degree lt k on triangle K P K all polynomials of total degree k on triangle K py e all polynomials of degree lt k on edge e Raviart Thomas space on triangles RT K P K 2P K k gt 0 dim RT k 1 k 3 First we can construct a basis for RT K i e we have RT K span w i 1 k 1 k 3 with J Uy 0 lt r s lt k ift 1 k 1 k 2 2 Bey 4 Ge yOsrts lt sk ifi K4 1 K4 2 241 6 1 k 2 3 eryar o lt r lt k if i k 1 k 2 41 k 1 k 3 4 For the construction of a H div Q conforming space we have to ensure continuity of the normal component across element edges This can be achieved by defining moments m and basis functions which are orthonormal with respect to this moments Let Ja nop ds pr pk eo O lt r lt k ifi 1 k 1 x Ja napi ds pr peer O lt r lt k ifi k 2 2k 2 mild J n p rds p pk e2 O lt r lt k ifi 2k 3 3k 3 fe 1 272. 13 13 14 14 14 14 15 15 15 15 8 1 Bilinear forms 8 2 Rieht hand sides r qq ll aha wh eR A Ie ee 9 Adaptive schemes 9 1 Storage for indicators and refinements 2000 4 9 2 Indicator components 10 Local basis functions 10 1 Scalar basis functions 2 2 a a mn e ee 10 2 Vector valued basis functions 2 0 0 ee ee 10 2 1 Raviart Thomas elements RT on squares 10 2 2 Raviart Thomas elements RT on triangles 10 2 3 Nedelec elements ND on cubes 2 222 ae 10 2 4 Nedelec elements ND on tetrahedrals 11 Polynomials 12 Interfaces to the integral libraries 12 1 2D BE Mince cs do eed ATA A RN a a a da ae 132 3D BEM u SE Ma Mat teh Dum Tae Nee a ey ee Ge h fan a gs 12 3 2D FEM au er 2020 0 2 ir ar Matam EA b 3 12 4 3D EEM 2 og aye ese se usata ete ea ee eine 13 Interface of the batch control language 14 Logging facilities 15 Parallelization 16 Internal organization of solvers 17 Internal organization of preconditioners 18 Configuration 19 Functions and boundary data 17 17 18 19 19 20 20 21 22 23 23 26 26 34 36 36 37 40 40 41 41 44 44 1 Introduction maiprogs is a fully fledged f90 program and makes use of most modern f90 features Espe cially the whole package is structured by the use of modules encapsulation by using generic interf
3. fx contains the result 44 of the evaluation of the function which is denoted by ltyp 0 given by the common block Itype There are printing functions for boundary data and volume data which describe the function in clear text on standard output subroutine pr op f ltyp integer intent in ltyp select case ltyp case 0 case 1 case default end select end subroutine pr op f subroutine prv op f ltyp integer intent in ltyp select case ltyp case 0 case 1 case default end select end subroutine prv op f References 1 M Maischak The analytical computation of the Galerkin elements for the Laplace Lam and Helmholtz equation in 2D BEM Preprint 95 15 DFG Schwerpunkt Randele mentmethoden 2 M Maischak The analytical computation of the Galerkin elements for the Laplace Lam and Helmholtz equation in 3D BEM Preprint 95 16 DFG Schwerpunkt Randele mentmethoden 3 M Maischak Manual of the program system maiprogs Preprint Institut fiir Ange wandte Mathematik Universitat Hannover 45
4. 0 lt r s lt k 1 i 3k 4 3k 3 k k 1 2 fe d ry0 lt r s lt k 1 ifi 3k 4 k k 1 2 k 1 k 3 5 Then the basis functions have to satisfy the condition mj fij We can represent i with respect to Di i e k 1 k 3 Y auth l 1 21 This leads to the following linear system k 1 k 3 mj i gt am D bi l 1 If we define bj m 4 then we have au 6305 For computational efficiency we later extend the coefficients ap with respect to 4 to a basis of P 1 K 2 which eliminates the need to distinguish several special cases in the evaluation of Qi 10 2 3 Nedelec elements ND on cubes Qim n all polynomials of degree lt in x direction of degree lt m in y direction and of degree lt n in z direction Nedelec space on cubes NDx Q u ur Qk 1 k k U2 Qk k 1 k U2 Qk k k 1 dim ND 3k k 1 First we have to construct a basis for ND Q i e we have ND Q span y i 1 3k k 1 2 with 1 O xryH 0 lt r lt k 1 0 lt s8s lt k 0 lt t lt k ifi 1 k k 1 0 _ 0 Yi x y z L atyszt O lt r lt k 0 lt s lt k 1 0 lt t lt k ifi k k 1 1 2k k 1 0 0 0 xrys2H 0O lt r lt k 0 lt s lt k 0 lt t lt k 1 ifi 2k k 1 1 3k k 1 1 6 For the construction of a H curl Q conforming space we have to ensure continuity of the tangential components across element faces This can be achieved by defining moments m and basis functions
5. boundary data traction Newton Potenial Poincare Steklov operator exterior domain Poincare Steklov operator interior domain inverse Poincare Steklov operator exterior domain Non linear term Matrix vector multiplication Spline will be tested 9 1 Storage for indicators and refinements The error indicators will be computed by the commands adap resh etc The error indicators itself will be stored in an real array with the name ekom in the bcdat database of the corre sponding spline Then the refinement information later used by refine will be computed and stored in the integer array with the name ref in the bcdat database of the same spline There can be different error indicators for independent splines which will be later refined independently 17 9 2 Indicator components Indicator Variable RES2CLAPUF u ETALAPUF Y hp llAun lar TETR 0 RES2CNEU u ETANEU mM 3 hellto Em an lirio Or RES2CJUMP u ETAJUMP Y he qe gt eedTn 1 1 RES2CTRANS u D N TUWKS n Y hellto Vun n 4 zW uo un ES Ten to IlZ2 e ui sa d IKVS m gt hell Z K uo un V n to Z RES2CDIVPF p DIVPF Mav gt divpn fler TET RES2CCURLP p CURLP Neat gt Ihr curl pallzacry TETn RES2CNORMP p NORMP m gt lihro TETR RES2CPTJUMP p PTJUMP Mo gt Whe pa t 22 o RES2CPTFREE p PTFREE m gt ht pn tlee ent pz RES2CGMPN u S N GMPN Mg log 1 Cz Em MP 9 pr li
6. exception handling incorporated into the language Then logmsg will evaluate the type argument for an extended action and error reporting The following subroutine timelog writes the consumed cpu time since the last command or timer reset to the screen and or log file The measured time will be marked by str If str is missing the current bcl command name will be used subroutine timelog str character len optional intent in str 15 Parallelization maiprogs currently supports three kinds of parallelization First parallelization on shared memory systems with OpenMP is fully supported and activated by default config openmp yes Second parallelization on clusters multi computer systems with MPI is supported partly still in development and has to be explicitly activated using config mpi yes Third parallelization using graphic cards with CUDA is in the early stages of support and has to be explicitly activated using config cuda yes 40 16 Internal organization of solvers The generic and reusable solvers are contained in the files fox lgsrev f90 for real numbers and fox xlgsrev f90 for complex numbers These solvers are used by the subroutines in the files fox lgs f90 which extend the solvers with the matrix vector multiplication for real dense matrices fox xlgs f90 with matrix vector multiplication for complex dense matrices fox femlgs f90 with matrix vector multiplication f
7. intent in ckom 0 x type spline1 intent in sp real kind dp intent out fwx 0 5 2 Two dimensional splines Compute the value of solution ckom at point x in element i 2D subroutine fpsi2 x i typ sp ckom fw al a2 a12 use poly only ptab2n basenum2 integer intent in typ i real kind dp intent in ckom 0 x 0 1 type spline2 intent in sp real kind dp dimension 0 sp bmode 1 1 intent in optional al a2 a12 real kind dp intent out fw 0 Compute the gradient of solution ckom at point x in element i 2D subroutine grdpsi2 x i typ sp ckom fwx fwy al a2 a12 use poly only pstab2n basenum2 integer intent in typ i real kind dp intent in ckom 0 x 0 1 type spline2 intent in sp real kind dp dimension 0 sp bmode 1 1 optional al a2 a12 real kind dp intent out fwx 0 fwy 0 5 3 Three dimensional splines Compute the value of solution ckom at point x in element i 3D subroutine fpsi3 x i typ sp ckom fw a 11 use poly only ptab3 basenum3 integer intent in typ i real kind dp intent in ckom 0 x 0 2 type spline3 intent in sp real kind dp intent in optional a 0 2 0 7 real kind dp intent out fw 0 Compute the gradient of solution ckom at point x in element i 3D subroutine grdpsi3 x i typ sp ckom fwx fwy fwz a use poly only pstab3 basenum3 integer intent in typ i real kind dp inten
8. solution with two or three components Note vtypc and vtypf can be different ng denotes the actual number of mesh elements ri i i 0 ng 1 gives the beginning of the local degrees of freedom in a global setting therefore ri i 1 ri i i 0 ng 1 gives the number of local degrees of freedom defined on the i th mesh element rci maps the local degrees of freedom to the global degrees of freedom i e do i 0 ng 1 do j 0 ri i 1 ri i 1 write i j x rci ri i j ktyptr r 0 ktyp 1 end do end do prints the values of the vector x for all local degrees of freedom rci is also used to eliminate degrees of freedom to incorporate homogeneous boundary conditions The corresponding local degrees of freedom are just mapped to the first element behind the global degrees of freedom The number of existing global degrees of freedom is obtained by dof In previous versions the global degrees of freedom could be obtained by n rci ri ng but this is not longer permissable 10 5 1 One dimensional splines Compute the value of solution ckom at point x in element i 1D subroutine fpsil x i sp ckom fw use poly only ptabi integer intent in i real kind dp intent in ckom 0 x type spline1 intent in sp real kind dp intent out fw 0 Compute the gradient of solution ckom at point x in element i 1D subroutine grdpsil x i sp ckom fwx use poly only pstabl integer intent in i real kind dp
9. Here vt 0 means the variable is of integer type with value vi and vt 1 means the variable is of real type with value vr If vt 1 the variable name is not defined subroutine getglobvar name vt vi vr integer intent out vt vi real kind dp intent out vr character len name The following subroutines and functions are designed for internal use ausset looks whether a file has been opened for output and if not it opens the file subroutine ausset ausgabe ausflag aus integer intent inout ausflag integer intent in aus character len intent in ausgabe iscmd is a logical function tests for the next command in the input line and compares with str The return value of iscmd is true after a command is recognized Then the string currentcmd will be set to the value of str and the subroutine timing will be called function iscmd str logical iscmd character len intent in str Example 13 1 This example shows the use of getopen getline and getclose I recin is the name of the default input file the main loop starts here note common batch commands are contained in the file bcl f and executed by the command getline call getopen recin recout fin 10 call getline fin if fin eq 1 then done don t do any more else if iscmd px then write mesh to output file else if iscmd mesh then initialisize the mesh call t
10. Technical Manual of the program system mazprogs Matthias Maischak Brunel University Uxbridge UK Institut fiir Angewandte Mathematik Universitat Hannover April 25 2015 O Overview This manual contains the parameter lists of subroutines in the libraries and other source files of the program system maiprogs which are shared by the different programs and which can also be used for own programs The interfaces are mostly stable and probably will not change in the future But its recommended that you try to use the newest version if possible For most internal data structures access functions are provided which will only enhanced in their functionality but which will remain compatible on the language level It is strongly recommended that you use data structures only via their access functions to keep the amount of effort limited which is necessary to update to a new maiprogs version For simple changes to maiprogs the config script allows to define additional directories in depmod conf which are used to search for main programs and additional modules and whose contents also override the versions in the main directories if you use the same name for a file like in the main directory By using this mechanism there is no need to change the original file at all Only make a copy to your own directory and change the copy Be aware you should never use the same name for a module in your own source file as in a standard source file when the sourc
11. aces and declaring the module procedures private and also by the definition of data types dedicated to the description of meshes and splines Global constants are always de clared in the module header Especially we have to note the definition of the kind parameter dp kind 0 0d0 which is used to define the precision of real variables which is part of the module const It is strongly recommended that you never access the internal data structures directly but use the access functions provided If for some reason you see no possibility to avoid the direct use of internal data structures either for reasons of speed or because some functionality is missing please contact the author of this package for providing the needed functionality in a future release 2 Components of a FEM BEM programme In this section we want to outline the general structure of a FEM BEM programme i e we want to present the sequence of steps which are to be taken to achieve an efficient FEM BEM code and give the reasons for the most crucial choices We do not want discuss here the general bookkeeping but we want to concentrate on the way the whole programme is broken down in smaller units which are mostly independent of each other and communicate only by some data structures There are the following components Geometry Mesh Splines Matrix Block Matrix Right hand Side Solver 3 Internal data structures for geometry definition There is only one general data fo
12. alue rvalue integer mode ivalue real rvalue where ivalue is an integer parameter and rvalue a real parameter mode 0 sets the kind of integration i e ivalue 0 means analytical integration ivalue 1 means doing the outer integration numerically and the inner integration analytically and ivalue 2 means all inte grations numerically The number of Gauss nodes gqna for the outer integration is set by mode 1 and the number of Gauss nodes gqnb for the inner integration by mode 2 The number of Gauss nodes gqnc for the potential integration is set by mode 3 sigma and ijn give the grading parameter and number of levels for the geometrical refined quadrature mu describes the increasing of the number of quadrature points farreg gives the boundary of the farfield where the usual quadrature rule is used without grading By using mode 1 the standard values are set this subroutine is automatically called by initgauss12 mode ivalue rvalue 0 mtyp 1 gqna 2 gqnb 3 gqne 4 ijn sigma 5 farreg 6 mu We can use the subroutine getc12 to obtain the current values of the parameters subroutine getc12 mode ivalue rvalue integer mode ivalue real rvalue In the following the interfaces and data formats of all user usable subroutines are presented The analogous subroutines of liblap2 f90 liblame2 f90 and libhelm2 f90 share a common interface Therefore we need only to state a generic interface The actual names of the subroutines are build by u
13. des px py pz typ ptyp real kind dp intent inout gm 0 real kind dp intent in nodes 0 2 0 7 integer intent in px py pz typ ptyp subroutine type 3 op gm nodes px1 py1 pz1 px2 py2 pz2 typ ptypl ptyp2 real kind dp intent inout gm 0 real kind dp intent in nodes 0 2 0 7 integer intent in px1 py1 pz1 px2 py2 pz2 typ ptyp1 ptyp2 36 13 Interface of the batch control language All programs of maiprogs read the commands from a file The implemented commands are described in the user manual 3 The subroutines and functions which are needed for implementing the interpretation of the commands are described in the following We have to distinguish between the routines which open and close the files and the routines for reading the parameters of commands The first subroutine is getopen which reads the arguments of the command line of the program itself o 02 f d i r and then opens the command file and also initializes most of the data structures subroutine getopen eingabe0 ausgabe0 fin only one call at program start for initialization integer intent inout fin error code emergency stop character len intent in eingabe0 ausgabe0 Here eingabe0 is the name of the default input file whereas ausgabe0 is the name of the default output file fin is a flag for emergency exit The second subroutine is getline which reads the next line from the input and tries to interpret the
14. dpoint diffvec and normal of boundary element 0 real m0 0 1 d0 0 1 n0 0 1 c function external fkt c maximal polynomial degree integer p0 c quadrature integer gqn real gqk 1 gqn gqw 1 gqn integer nt opt c singular point real sx 0 1 c table of the identity integer fm 33 parameter fm 20 real kid 0 fm 0 nt 1 12 2 3D BEM The source code files of the integral libraries in the 3D case are called liblap3 f90 liblame3 f90 libhelm3 f90 libint3 f90 and libgint3 f90 The file liblap3 f90 contains the implementation of the Galerkin elements for the Laplacian the file liblame3 f90 contains the implementation of the Galerkin elements for the Lam operator and the file libhelm3 f90 contains the im plementation of the Galerkin elements for the Helmholtz operator Each of this files can be used independently of the others but all files need the file libint3 f90 which contains the elementary integrals and the file libgint3 f90 which contains generic numerical integration subroutines used by the other files Before using any subroutines in which numerical integrations are involved we have to call the subroutine initgauss23 with a subroutine as argument which calculates for a given number of nodes the values of the nodes and weights The subroutines initgauss23 setc23 and getc23 correspond to the subroutines initgauss12 setc12 and getc12 in the 2D BEM case and share the same interface subroutine integ3 bx dx1 dx2 nx tx by dy1 dy2 n
15. e files are differently named This will confuse the build mechanism Remark There is no warranty that these routines will work correctly as specified The use of these routines is on your own risk Using these routines in commercial programs is prohibited as long as not explicitly allowed by the author email matthias maischak brunel ac uk Contents O Overview 1 Introduction 2 Components of a FEM BEM programme 3 Internal data structures for geometry definition 3 1 Geometry database ee ee ee 4 Internal data structures for Meshes 4 1 One dimensional manifolds 2 22 22 2 no nn e e 4 2 Two dimensional manifolds nn nn 4 3 Three dimensional spaces e 4 4 Mesh databases 0 nen 5 Internal data structures for Splines 5 1 One dimensional splines 2 00000 eee ee ee 5 2 Two dimensional splines 2 22 22 2 CE mE 5 3 Three dimensional splines 2 2 2 2 on on nn 5 4 Spline databases naa br ki doe ek er aa 5 9 bedat database san ana b Yay ur a Er kas s es a ni 6 Problem data structure Gall o SS An E eek ee 6 2 datstry ara A A re a Ale ate gen ee 6 3 SOOSTE 2 23 a ee ed hates 7 Matrices TT Sparse Matrices usada Sah bd ee Be ee GG 1 27 Dense Matrices tar tot Bin Bee Fa Be IS 3 E M trices and fake See et pO qaq en Sa EG TA Block Matrices aqa s ae o A a E Pe RES 8 Bilinear forms and right hand sides 11 11 11 12 13
16. e sep arated by or All matrices after a are still in the matrix database and will be managed and computed automatically but are not part of the system matrix The entries have the format test spline Matrix name arguments trial spline e g u A v t u for a ma trix depending on the splines v and t where test and trial splines are the same matstr 13 is analyzed by the subroutine setmatrix the Matrix name is evaluated by the subroutine matcomp2c matcomp3c etc With defmatrix we assign a bilinear form to the matrix A list of available bilinear forms is given in Section 8 1 6 2 datstr Using datstr we describe the construction of the right hand side Entries are separated by The entries have the form test spline expression expression can be any linear com bination of given functions f u0 t0 given splines operators applied to functions or splines and matrices multiplied with splines datstr is analyzed by the subroutine setlft A list of available right hand side operations is given in Section 8 2 6 3 genstr The order in genstr determines the order in which the coefficients are ordered in the coefficient 9 vectors ckom and lkom Entries are separated by or Only one is allowed Entries before belong to the unknowns in a linear equation system whereas entries after are used to denote auxilliary splines The following fu
17. eger pl integer nt opt c table of the integrated potential integer fm parameter fm 20 real potkr 0 fm 0 nt 1 Numerical quadrature of a numerical potential with a function subroutine genpotf1n m0 d0 n0 m1 d1 n1 f kernf mult pi potkr nt opt c midpoint diffvec and normal of boundary element 0 real m0 0 1 d0 0 1 n0 0 1 c midpoint diffvec and normal of boundary element 1 real m1 0 1 d1 0 1 n1 0 1 c external subroutine with parameters x n f 32 c x 0 1 source point c n 0 1 normal direction in x c 0 1 values in x external f kernf mult c maximal polynomial degree integer pl integer nt opt c table of the integrated potential integer fm parameter fm 20 real potkr 0 fm 0 nt 1 Numerical quadrature of a double integral consisting of a singular kernel and two functions subroutine gengalf1n m0 d0 n0 m1 d1 n1 f0 f1 kernf emult galr nt opt c midpoint diffvec and normal of boundary element 0 real m0 0 1 d0 0 1 n0 0 1 c midpoint diffvec and normal of boundary element 1 real m1 0 1 d1 0 1 n1 0 1 external subroutine with parameters x n f x 0 1 source point n 0 1 normal direction in x f 0 1 values in x external f0 f1 kernf emult c maximal polynomial degree integer pl integer nt opt 000 0 c table of the integral real galr 0 nt 1 Computation of the identity where f is allowed to have a singular point hen f ro Saca subroutine genid12 m0 d0 n0 p0 gqn gqk gqw kid fkt nt opt sx c mi
18. ement 1 real m1 0 1 d1 0 1 n1 0 1 integer nt c kernel function arbitray function elementary multiplication external kernf fkt emult c table of the potential values real potw O nt 1 Numerical computation of the Galerkin elements given only the kernel function subroutine gengalln m0 d0 n0 m1 d1 n1 p0 p1 f kernf nt opt c midpoint diffvec and normal of boundary element 0 31 real m0 0 1 d0 0 1 n0 0 1 c midpoint diffvec and normal of boundary element 1 real m1 0 1 d1 0 1 n1 0 1 integer p0 p1 nt opt integer fm parameter fm 20 real f 0 fm 0 fm 0 nt 1 Numerical computation of the outer integration given the analytical computed potential subroutine gengall m0 d0 n0 m1 d1 n1 p0 p1 f potl11 nt opt c midpoint diffvec and normal of boundary element 0 real m0 0 1 d0 0 1 n0 0 1 c midpoint diffvec and normal of boundary element 1 real m1 0 1 d1 0 1 n1 0 1 integer pO pi integer fm parameter fm 20 real f 0 fm 0 fm 0 nt 1 external potll Numerical quadrature of a potential with a function subroutine genpot m0 d0 n0 m1 di ni f potll mult pi potkr nt opt c midpoint diffvec and normal of boundary element 0 real m0 0 1 d0 0 1 n0 0 1 c midpoint diffvec and normal of boundary element 1 real m1 0 1 d1 0 1 n1 0 1 external subroutine with parameters x n f x 0 1 source point n 0 1 normal direction in x 0 1 value in x external f potll mult c maximal polynomial degree Qo Q QO O int
19. er ng bmode 0 1 character len 32 name end type meshi The mesh consists of ng boundary elements given by start points rx i and distance vectors rh i i 0 ng 1 For every boundary element the normal direction is given by rn i rc 0 i gives the element number of the left neighbor connected to the point rx i and rc 2 i gives the element number of the right neighbor connected to the point rx i rx i In most cases we don t need to deal with the internals of mesh representation The access subroutine subroutine mh2mdn mh i mx dx nx type mesh1 intent in mh integer intent in i real kind dp intent out mx 0 1 dx 0 1 nx 0 1 extracts the midpoint mx 0 1 the distance vector dx 0 1 pointing from the midpoint to the right endpoint and the normal direction nx 0 1 of the i th element from the mesh mh independent of the internal representation 4 2 Two dimensional manifolds type elem2 integer rnodes 0 4 rnodes 0 3 triangle 4 quadrilateral rnodes 1 4 node numbers integer rc 0 3 neighbours integer rct 0 3 sides of neighbours integer father integer sons 1 4 integer ref real kind dp rn 0 2 Normal for 3d mesh end type elem2 type mesh2 real kind dp dimension pointer rx rdi rd2 rn integer dimension pointer rt integer dimension pointer rc rct real kind dp dimension pointer nodes integer nnode
20. ext to standard output by printpoly where fin is an error code if a wrong non existent type is chosen subroutine printpoly ptyp fin integer ptyp fin Other subroutines recommended for public use are the following Print the internal one dimensional polynomial coefficients of order i to channel out subroutine printpmi i out use polydat only npmi pmi integer intent in i out Evaluate base functions up to degree p at point x subroutine ptabl x p ptyp tab integer intent in p ptyp real kind dp intent in x real kind dp intent out tab 0 max 1 p Evaluate derivatives of base functions up to degree p at point x subroutine pstab1 x p ptyp tab use polydat only pmi npmi integer intent in p ptyp real kind dp intent in x real kind dp intent out tab 0 max 1 p 23 Evaluate second derivatives of base functions up to degree p at point x subroutine psstab1 x p ptyp tab use polydat only npmi pmi integer intent in p ptyp real kind dp intent in x real kind dp intent out tab 0 max 1 p Evaluate base functions up to degree px py at point x y subroutine ptab2n x y px py typ ptyp tab use polydat only tslegs real kind dp intent in x y integer intent in px py typ ptyp real kind dp intent out tab 0 Evaluate derivatives of base functions up to degree px py at point x y tabl contains the x derivative and tab2 the y derivative sub
21. f90 libhelm2 f90 libint2 f90 and libgint2 f90 The file liblap2 f90 contains the implementation of the Galerkin elements for the Laplacian the file liblame2 f90 contains the implementation of the Galerkin elements for the Lam operator and the file libhelm2 f90 contains the implementation of the Galerkin elements for the Helmholtz operator Each of these files can be used independently of the others but all files need the file libint2 f90 which contains the elementary integrals and the file libgint2 f90 which contains the generic numerical integration subroutines used by the other files Before using any subroutines in which numerical integrations involved we have to call the subroutine initgauss12 with a subroutine as argument which calculates for a given number of nodes the values of nodes and weights initgauss12 will generate internally a table for the number of nodes from 1 to 32 for further use subroutine initgauss12 gauss external gauss An example for such a gauss routine is given in the following by using NAG routines subroutine gauss gqn gqk gqw integer gqn itype ifail real gqk 1 gqn gqw 1 gqn real a b external d0ibaz d01bbf a 1 0 b 1 0 itype 0 ifail 1 call dO1bbf d01baz a b itype gqn gqw gqk ifail if ifail ne O then 26 ifail 0 call dO1bcf itype a b a b gqn gqw gqk ifail end if end Internal coefficients controlling the analytical or numerical integrations are set by subroutine setc12 mode iv
22. i y 2 lt l lt p left edge bri y Lrr la y 1 lt k lt p 1 right edge oril y Leoa a i y 1 2 y 1 lt k 2 lt 1 k 1 lt p_ inner bubble functions For the construction of continuous splines we have to take into account the orientation of the edge functions ie 1 10 2 Vector valued basis functions 10 2 1 Raviart Thomas elements RT on squares Qim all polynomials of degree lt l in x direction and of degree lt m in y direction p e all polynomials of degree lt k on edge e Raviart Thomas space on squares RTy K Qk 1 k x Orjat dim RT K 2 k 1 k 2 First we can construct a basis for RT K i e we have RT K span w i 1 2 k 1 k 2 with gt JS JOYO0SEr lt k 10 lt s lt k ifi 1 k 1 2 Way S atys O lt r lt k O lt s lt k4 1 ifi k 1D k 2 1 208 1 k 2 2 For the construction of a H div Q conforming space we have to ensure continuity of the normal component across element edges This can be achieved by defining moments m and basis functions which are orthonormal with respect to this moments Let Sts ms lt ifi k 2 2k 2 e2 O lt r lt k ifi 2k 3 3k 3 B 3 Mm 3 m wo o A 3 A ifi 3k 4 4k 4 ifi 4k 5 k 1 k 4 ifi k4 1 k 4 1 2 k 1 k 2 3 20 Then the basis functions i have to satisfy the condition mj i fij We can represent i with respect to Dr i e
23. imelog write computing time to screen and or logfile else error message call errcmd end if if fin ne 1 goto 10 39 call getclose 14 Logging facilities Messages informational warnings errors are written to the screen and or to a log file In the module error there is the subroutine logmsg and the character variable msg containing the message subroutine logmsg type name only character len intent in optional type name only Supplying the optional argument type we can define the type of the message e g Error Warning Info Supplying the optional argument name we can use the name of the reporting procedure or command in the message This has the advantage that the reporting procedure can call first another subroutine for further processing supplying the name which finally leads to the call to logmsg Using the optional argument we can supercede the general configuration writing only to Screen or File It is strongly advised to use the logging facilities for all kinds of structured output i e for anything besides debugging outputs so that the output can also be organized in a multi core multi computer environment We have to note that the integral libraries are not covered by the logging facilities The integral libraries are intended to be used also externally therefore a close connection would be harmful for this purpose In the next Fortran revision there will also be some kind of
24. line as a standard command in out out2 e set print write2 do continue history debug id ti time If it recognizes no command the input line is returned to the calling program for further interpreting It also controls the history buffer do loop control etc subroutine getline fin integer intent inout fin The last subroutine is getclose which just closes input and output files subroutine getclose Now the subroutines and functions for interpretation of the input lines are given All com mands which read from the input line use the pointer inpp internally outset looks whether a channel nr 1 upto 10 has been opened for output and if not it opens the corresponding file subroutine outset chno aus integer intent in akno integer intent out aus create a new file for channel akno allocated unit will be in aus rdint reads an integer from the input line an argument has to be present otherwise an error message will be given function rdint integer rdint rdreal reads a real number from the input line an argument has to be present otherwise an error message will be given 37 function rdreal real kind dp rdreal xrdint reads an integer from the input line and uses a default value wert if there is no more input or the symbol is used function xrdint wert integer xrdint integer intent in wert xrdreal reads a real number from the input line and uses a defau
25. lli do 2 for all k l do aster fa Vv ee end for for all k do p lt Tik for alll do q Til Ap Qa 2 wi pwi ae Em 10 end for 11 end for 15 12 end for 8 1 Bilinear forms We have to distinguish bilinear forms which lead to sparse matrices from bilinear forms which lead to dense matrices They differ in the way they are assembled Laplace We have the following implemented Bilinear forms for FEM and BEM A M B AO Al HDIV PenaltyO DGE DGF DGC V K K LK LI K Qij Vo dx Qij 6 0 dx Q Qij div p u dx Q Q Q Single layer potential Vy 2 ete Gle most ds ds Double layer des Ky 2 dle E G x y Y u ds ds W W PP Hypersingular creel operator S Exterior Poincare Steklov Operator S W K D V K D Sin Interior Poincare Steklov Operator S m K NVUK N R Inverse Exterior Poincare Steklov Operator Lame Qij 000 dx Helmholtz A Qij Vo Vo dx Bid dx Q Q M ay 6 6 dz Q Maxwell Convection CONV Galerkin Matrix SUPG SUPG Matrix CMDir 16 Stokes M VDiv STr V V NN K K LK LI K W W PP NN NN bier 07 1 Bilaplace A M ay oud dx 8 2 Right hand sides u0 t0 NV NK SUPG VW K Ks K I K Ks 1 Ks I S Sin R NLin Matrix Spline Spline 9 Adaptive schemes m m bier ds moja dsy x di ass nulo dsy Volume data boundary data displacement
26. lt value wert if there is no more input or the symbol is used function xrdreal wert real kind dp xrdreal real kind dp intent in wert The most modern and comfortable way to access arguments for bcl commands is to use the following set of commands yparse yrdint yrdreal and yrdstr The subroutine yparse reads all arguments at once This makes it possible to enclose the arguments in brackets and to separate them with colons which was not possible before Additionally the arguments can be named and therefore given in arbitrary order subroutine yparse The name of the optional argument is given in str and its default value in value function yrdint str value integer yrdint integer intent in value character len intent in str function yrdreal str value real kind dp yrdreal real kind dp intent in value character len intent in str function yrdstr str value character len 80 yrdstr character len intent in str value The subroutine setglobvar allows to set global variables which can be used by bcl scripts Here vt 0 means the variable is of integer type with value vi and vt 1 means the variable is of real type with value vr subroutine setglobvar name vt vi vr integer intent in vt vi real kind dp intent in vr character len intent in name 38 The subroutine getglobvar allows to get global variables which are used by bcl scripts
27. me files The three types splinel spline2 and spline3 have the following components in common e The component ng denotes the number of elements of the mesh to which the spline is subordinated e The component name contains the full name of the spline variable e The component bmode 0 1 contains in bmode 0 the dimension of the manifold and in bmode 1 the dimension of the space Initializing the internal coefficients necessary for a spline defined in the module polydat can be done by calling the generic subroutine subroutine polyinit sp use spline 1 2 3 type spline 1 2 3 intent in sp The individual representations are defined in splinel spline2 and spline3 This section deals with the most important data structures of the whole program system i e the structures for meshes and solutions We use a consistent scheme for all cases 2D BEM 3D BEM 2D FEM but note that our scheme is most appropriate for the implemen tation of a general hp mesh i e an arbitrary distribution of mesh elements and polynomial degrees Lets first make a few restricting assumptions on the kind of splines which we want to deal with Later we can generalize everything Our fundamental assumptions are the following e Our global basis functions are based on a mesh e On every mesh element we have a given set of local basis functions e The local basis functions are generated by mapping a reference element to the mesh element e Basis func
28. n x and nl in y kern O nt 1 is the value of the kernel function with nt components For the implementation of integral equations of the second kind we have to test a base function with another base function Ia Lor 91 1 yond N e 7 de des To Ty This is given by the subroutine genidgal 30 subroutine genidgal m0 d0 n0 m1 d1 n1 p0 p1 idk1 c midpoint diffvec and normal of boundary element 0 real m0 0 1 d0 0 1 n0 0 1 c midpoint diffvec and normal of boundary element 1 real m1 0 1 d1 0 1 n1 0 1 integer pO pi integer fm parameter fm 20 real idk1 0 fm 0 fm Various subroutines located in libgint2 f90 perform the numerical quadratures Numerical computation of the potential for monomials given only the kernel function subroutine genpotln x n0 m1 d1 n1 py potkr kernf nt opt c source point of the potential real x 0 1 c normal of boundary element 0 real n0 0 1 c midpoint diffvec and normal of boundary element 1 real m1 0 1 d1 0 1 n1 0 1 c maximal polynomial degree integer py nt external kernf c tabler of the values of the potential integer fm parameter fm 20 real potkr 0 fm 0 nt 1 Numerical computation of the potential for arbitrary functions given only the kernel function subroutine genfpotln x n0 m1 d1 n1 potw kernf fkt emult nt opt c source point of the potential real x 0 1 c normal of boundary element 0 real n0 0 1 c midpoint diffvec and normal of boundary el
29. n0 al el n1 f n op kr c startpoint endpoint and normal of the O boundary element real a0 0 1 e0 0 1 n0 0 1 c startpoint endpoint and normal of the 1 boundary element real a1 0 1 e1 0 1 n1 0 1 29 external subroutine with parameters x n f x 0 1 source point n 0 1 normal direction of the boundary in x 0 1 value of the function in x external f c highest polynomial degree integer n oo y O c table of the values integer fm parameter fm 20 real op kr 0 fm 0 nt 1 These subroutines also exist using midpoint and difference vector subroutine op potmd m0 d0 n0 m1 d1 n1 f n op kr The values of the potential at a given point Fa FEDD As Ti can be computed by using the following subroutines which use the same parameter list subroutine op pot11 x n0 m1 d1 n1 p1 op kr c source point real x 0 1 c midpoint diff vector and normal direction of the boundary element 1 real m1 0 1 d1 0 1 n1 0 1 c normal direction in the source point x real n0 0 1 c highest polynomial degree integer pl c table of the values integer fm parameter fm 20 real op kr 0 fm 0 nt 1 The value of the potential can be calculated by summing up the results of op potll over the whole boundary The kernel functions F x y of the integral operators are given by subroutine op kern d n0 n1 kern real d 0 1 n0 0 1 n1 0 1 kern 0 nt 1 d 0 1 is the difference of x and y n0 denotes the normal direction of the boundary i
30. nctionality is provided here u p denote splines in the domain and D D1 splines on the boundary but the names are fully arbitrary but have to be chosen globally unique u The Spline space u will be initalized and the memory for a Spline called u will be allocated u p The Spline space u will be initalized based on the mesh used by p and using the same polynomial degrees as p u p 1 The Spline space u will be initalized based on the mesh used by p and using the polynomial degrees of p plus 1 u 1 The Spline space u will be initalized using the previously given mesh and the polynomial degree 1 u Spline p The memory for a spline called u will be allocated using the Spline space p D Dir u D will be the trace space of u the coefficients of u will be reordered such that D will have consecutive coefficients D lt Dir u 1 D will use the trace mesh of u and polynomial degrees plus 1 but will not share the coefficients with u Di lt Restrict D2 G1 1 D will use the mesh and polynomial degrees of D2 restricted to the geometry G1 but will not share the coefficients with D2 Di Reord D D1g Di lt Reord D D1g 1 7 Matrices 7 1 Sparse Matrices Sparse matrices are represented by the data type matrix_sparse The search function to access the sparse matrix database is qsmat i e with smat gt qsmat st
31. ngle the fourth node is an additional node on the surface The next three nodes define the corners of the second curved triangle the last node is an additional node on the second surface 3 1 Geometry database The search function to access the geometry database is qgm i e with gm gt qgm str the pointer gm contains the result of our database query for a geometry with the name given by the variable str Don t confuse the database of defined geometries with the database of pre defined geometries ads geometries 4 Internal data structures for Meshes There are three data types created for the description of meshes in defined on a manifold with 1 2 or 3 dimensions Note that the dimension of the manifold can be different from the dimension of the space Depending on this dimension the definitions are contained in the modules splinel spline2 or spline3 located in the file splines1 f90 splines2 f90 or splines3 f90 mesh and spline definitions are located in the same files The three types mesh1 mesh2 and mesh3 have the following components in common e The component ng denotes the number of elements of the mesh e The component name contains the full name of the mesh e The component bmode 0 1 contains in bmode 0 the dimension of the manifold and in bmode 1 the dimension of the space 4 1 One dimensional manifolds type meshi real kind dp dimension pointer rx rh rn integer dimension pointer rc rct integ
32. ntial which are usually computed by taking the transposed matrix of the double layer potential in the same way subroutine integmdks m0 d0 n0 m1 d1 n1 pks0 pks1 fks midpoint diffvec and normal of boundary element 0 real kind dp intent in m0 0 1 d0 0 1 n0 0 1 midpoint diffvec and normal of boundary element 1 real kind dp intent in m1 0 1 d1 0 1 n1 0 1 integer intent in pksO pksi integer parameter fm 40 real kind dp intent inout fks 0 fm 0 fm 0 nt 1 For the computation of the right hand side we need the integrals over a general function multiplied with the test functions Lx F Z Po 1 ds To subroutine id a0 e0 n0 f n gqn gqk gqw iid c startpoint endpoint and normal of the O th boundary element real a0 0 1 e0 0 1 n0 0 1 external subroutine with parameters x n f x 0 1 source point n 0 1 normal direction of the boundary in x f 0 1 value of the function in x external f GO O O c maximal polynomial degree integer n c quadrature integer gqn real gqk 1 gqn gqw 1 gqn c table of the identity integer fm parameter fm 20 real iid O fm 0 nt 1 For the computation of the right hand side we also need F l TG J PE Merit dsy dse which can be computed by using the following subroutines which use the same parameter list where op stands for v weakly singular k double layer ks adjoint double layer and w hypersingular e g lamevpot subroutine op pot a0 e0
33. op and the BPX algorithm a subroutine bpx mtop mdc xfm xfn x f pro prow difsi difsw drh rp ri ngm ngn nt ktyp computes one bpx step at level mtop integer nt ktyp mdc real f 0 x 0 real drh 0 2 0 integer rp O nt 1 0 ri 0 1 0 Position und Anzahl der Netze integer ngm 0 mtop ngn 0 mtop subroutine hiera ktyp mtop xfm xfn x f hpro hprow hdiag Projeziere auf die hierarchische Basis Teile durch die Diagonale Projeziere zurueck integer mtop ktyp integer xfm 0 mtop xfn 0 mtop real 0 x 0 hdiag 0 integer prom parameter prom 9 integer hpro 0 prom 0 real hprow 1 prom 0 43 18 Configuration makefile is the topmost Makefile whereas every folder contains its own Makefile which is called by the topmost Makefile config is a shell script which determines the operating system the used FORTRAN compiler and which libraries are available Using this infor mations config sets the Compiler options and link paths which are necessary to compile the whole program system It decides also whether it is necessary or possible to add some code to replace the BLAS and LAPACK libraries which are primarily used by an optimized version e g Intel s MKL this code is located in the directory ads It generates the files makefile mai and makefile gen 19 Functions and boundary data Functions describing volume data boundary data nonlinear behavior and exact solution
34. or real and complex sparse matrices and fox couplgs f90 with dense and sparse matrices for coupled systems 17 Internal organization of preconditioners Most preconditioners on matrix level are contained in the file fox mlbas f90 for dense sys tems and in fox mlfbas f90 for sparse systems Allimplemented preconditioners are of Schwarz type These preconditioners are distinguished by the construction of their subspace decomposition In general we can write that a space V with base is decomposed in subspaces V W V Vw Vi spanfd i 1 dimV For the implementation we have to describe the construction of amp in terms of j 1 The simplest subspace decomposition is given by a partition of the base functions Restriction and prolongation operators are simply given by the identity acting on a subspace 2 Another possibility is a subspace decomposition given by a partition of the base func tions with an additional coarse grid on which base functions are given as linear com binations of an arbitrary number of fine grid base functions These are the so called Hh methods 3 A famous possibility is nested subspaces where we have a sequence of even coarser subspaces grids which are constructed by the linear combination of few fine grid base functions to a coarse grid base function In the multiplicative Schwarz case these are the multigrid methods and in the additive Schwarz case the BPX like methods 4 A furthe
35. oundaries normal direction real kind dp dimension 0 2 intent in by dy1 dy2 ny Element typ 3 triangle 4 rectangle integer intent in ty 35 polynomial degree integer intent in py 0 1 table of the potential real kind dp intent inout op kl1 0 fm 0 fm 0 The value of the potential can be calculated by summing up the results of op potll over the whole boundary The kernel functions F x y of the integral operators are given by subroutine op kern d nO ni kern real kind dp intent in d 0 2 n0 0 2 n1 0 2 real kind dp intent out kern 0 nt 1 d 0 2 is the difference of x and y n0 denotes the normal direction of the boundary in x and nl in y kern O nt 1 is the value of the kernel function with nt components 123 2D FEM In case of the 2D FEM Galerkin elements we have to distinguish between Galerkin matrices with test and ansatz spaces which are of the same kind subroutine type 2 op gm p1 p2 p3 p4 px py typ ptyp use poly real kind dp intent inout gm 0 real kind dp dimension 0 1 intent in p1 p2 p3 p4 integer intent in px py typ ptyp and with different test and ansatz spaces subroutine type 2 op gm p1 p2 p3 p4 px1 py1 px2 py2 typ ptyp1 ptyp2 real kind dp intent inout gm 0 real kind dp dimension 0 1 intent in p1 p2 p3 p4 integer intent in px1 py1 px2 py2 typ ptyp1 ptyp2 12 4 3D FEM subroutine type 3 op gm no
36. pol integer intent in n i real kind dp intent in xi 0 n real kind dp intent out lagpol 0 n Compute uniform nodes for continuous ctyp 1 or discontinuous ctyp 0 interpolation on an interval a b subroutine ipnodes n ctyp a b xi integer intent in n ctyp real kind dp intent in a b real kind dp intent out xi 0 n Transformation from monomial base to base ptyp on an interval subroutine ctpol1 coeff p ptyp integer intent in p ptyp real kind dp intent inout coeff 0 max 1 p Transformation from monomial base to base ptyp on a rectangle or a triangle subroutine ctpol2 coeff px typ ptyp integer intent in px 0 1 typ ptyp integer parameter fm 10 real kind dp intent inout coeff 0 fm 0 fm Transformation of a polynomial with monomial base functions to a polynomial with Legendre polynomials as base functions subroutine ctlegi coeff p use polydat integer intent in p real kind dp intent inout coeff 0 p Transformation of Legendre coefficients to Antiderivatives of Legendre polynomials subroutine cbfleg1 coeff p integer intent in p real kind dp intent inout coeff 0 p 25 j 0 jz p L a n G z p p La n e a j 0 k 0 FAL ith JA p p p Kae gt SS ea 12 Interfaces to the integral libraries 12 1 2D BEM The source code files of the integral libraries are called liblap2 f90 liblame2
37. r bmode 0 1 integer nnodes nrnodeg ngg integer nrnodes ngs character len 64 name end type geom bmode 1 1 type nr 0 2 line segment bmode 1 2 type nr 0 2 line segment 1 3 arc segment described by three nodes located at the beginning the middle and the end of the segment 10 n Curvi linear element represented by Lagrange based uniform interpolation polynomial with n nodes bmode 2 2 type nr 0 3 triangle 0 4 quadrilateral element 1 4 triangle with one curved arc boundary segment the first three nodes define the arc segment whereas the last node defines a corner 1 5 quadrilateral element with one curved arc boundary segment the first three nodes define the arc segment whereas the two nodes define corners Effectively nodes 1 3 4 5 define the corners of the curved quadrilateral element whereas node 2 is used to define the curvature 1 6 quadrilateral element with two curved arc boundary segments opposite to each other Nodes 1 3 4 6 define the corners of the curved quadrilateral element Node 2 defines the curvature of the first arc segment and node 5 defines the curvature of the second arc segment 10 n Curvi linear element square represented by Lagrange based uniform interpolation polynomial with n nodes bmode 2 3 type nr 0 3 triangle 0 4 quadrilateral element 1 4 triangle with one curved arc boundary segment the first three nodes define
38. r the pointer smat contains the result of our database query for a sparse matrix with the name given by the variable str 14 7 2 Dense Matrices Dense matrices are represented by the data type matrix_dense The search function to access the dense matrix database is qdmat i e with dmat gt qdmat str the pointer dmat contains the result of our database query for a dense matrix with the name given by the variable str 7 3 H Matrices H matrices are represented by the data type matrix_hmat The search function to access the H matrix database is qhmat i e with hmat gt qhmat str the pointer hmat contains the result of our database query for a H matrix with the name given by the variable str 7 4 Block Matrices There exists a container datatype matrix_block which collects all the individual matrix formats Using bmat qblock str we search all the individual databases and obtain a container variable with pointers to any matrix found Note that bmat is not a pointer in itself 8 Bilinear forms and right hand sides The choice of global basis functions 1 has far reaching consequences Right hand sides of the form b Js f r p 0 dof 1 can be assembled in the following way 1 for all i do 2 for all k do 3 Wm fi JOSE 4 end for 5 for all k do 6 p Tik 7 b lt bp wide 8 end for 9 end for Sparse matrices of the form apq i VopVoq Pp q 0 dof 1 can be assembled in the following way 1 for a
39. r case is a complete transformation of the standard base functions to a base with special properties where usually a one dimensional Schwarz preconditioner the diagonal scaling is used Known examples are the hierarchical base and prewavelet or wavelet bases In the following we describe the interface to the implementations of these preconditioners 1 The partition of the base functions is given by the arrays integer nxm ixm 0 xm 0 nxm am 0 nxm An additional array for the local inverse is given by integer ipvt 0 Using these definitions we implement the additive Schwarz preconditioner 41 subroutine asmp ktyp b ixm xm am nxm ipvt n z r integer ktyp n real b 0 x z 0 n 1 r 0 n 1 the multiplicative Schwarz preconditioner subroutine msmp ktyp b a ix ixm xm am nxm ipvt flag n z r integer ktyp flag n ix 0 real a 0 x b 0 x z 0 n 1 r 0 n 1 and the symmetric multiplicative Schwarz preconditioner subroutine msmp2 ktyp b a ix ixm xm am nxm ipvt flag n z r integer ktyp flag n ix 0 real a 0 b 0 z O n 1 r 0 n 1 The coarse grid base is given by the arrays integer ncxm icxm 0 cxm O ncxm real wcxm 0 x An additional array for the coarse grid inverse is given by integer cipvt 0 Now the additive Schwarz preconditioner with coarse grid is given by subroutine varasmp ktyp b ixm xm am nxm ipvt icxm wcxm cxm ncxm cipvt n z r integer ktyp n real b 0 z 0 n 1 r 0 n 1 The memory di
40. red in the same data structure integer parameter fm 40 real kind dp 0 fm 0 fm 0 nt 1 fm denotes the highest allowed polynomial degree which can be used nt stands for the number of components of the Galerkin element For example for the Lam equation we have four components i e nt 4 We can access the Galerkin element F ie r s 0 1 by F k 1 2xr s8 The Galerkin elements are computed by integmd which uses midpoint difference vector and normal vector for the boundary elements 0 and 1 The kind of integration to be performed is defined in advance by setc12 0 ivalue rvalue The highest polynomial degrees on boundary elements 0 and 1 must be given for the single layer potential by pv0 and pvl for the double layer potential by pkO and pk1 and for the hypersingular operator by pw0 and pwl Note that not the hypersingular operator itself is computed but a partial integrated form which can be used completely analogous subroutine integmd m0 d0 n0 m1 d1 n1 pv0 pv1 pk0 pk1 pw0 pw1 fv fk fw midpoint diffvec and normal of boundary element 0 real kind dp intent in m0 0 1 d0 0 1 n0 0 1 midpoint diffvec and normal of boundary element 1 real kind dp m1 0 1 d1 0 1 n1 0 1 integer intent in pv0 pv1 pk0 pk1 pw0 pwl integer parameter fm 40 real kind dp dimension 0 fm 0 fm 0 nt 1 intent inout fv fk fw 28 We can compute the Galerkin elements for the adjoint double layer pote
41. riangle 4 rectangle integer intent in tx vector of the origin boundaries normal direction real kind dp dimension 0 2 intent in bx dx1 dx2 nx polynomial degree integer intent in px 0 1 quadrature integer intent in gqn real kind dp intent in gqk 1 gqn gqw 1 gqn integer intent in nt table of the identity real kind dp intent inout kid 0 fm 0 fm 0 nt 1 For the computation of the right hand side we also need Pex Ho N FDP de ds To Ty which can be computed by using the following subroutines which use the same parameter list where op stands for v weakly singular k double layer ks adjoint double layer and w hypersingular e g lamevpot subroutine op pot bx dx1 dx2 tx nx by dy1 dy2 ty ny f py op kr Element typ 3 triangle 4 rectangle integer intent in tx ty vector of the origin boundaries normal direction real kind dp dimension 0 2 intent in bx dx1 dx2 nx real kind dp dimension 0 2 intent in by dy1 dy2 ny polynomial degree integer intent in py 0 1 table of the potential real kind dp intent inout op kr 0 fm 0 fm 0 The values of the potential at a given point Fula f FED ds Ty can be computed by using the following subroutines which use the same parameter list subroutine op pot11 x nx by dy1 dy2 ny ty py op k1 Isource point real kind dp intent in x 0 2 nx 0 2 vector of the origin b
42. rmat for geometry definitions in maiprogs This general container format is used in all 3 dimensions The components ending on g contain a geometry definition which plainly defines the ge ometry itself The components ending on s contain an alternative definition of the same geometry with additional information on the location of singularities which is used for the generation of structured meshes In the following we will omitt the trailing g or s ng contains the number of macro elements in the geometry definition rnode the node numbers for each macro element type the kind of macro element rn an additional normal direction dn is a descriptor for every macro element which can take any value currently alway 0 ori in case of singularities describes the location node contains the coordinates of all nodes bmode 0 1 as usual describes the dimension of the geometry The one dimensional rnode structure can only be read consecutively For every element the first entry is the number of nodes nr belonging to this element followed by the node numbers type and number of nodes nr for every element define the kind of macro element type geom integer dimension pointer rnodeg typeg dng real kind dp dimension pointer rng integer dimension pointer rnodes types ori dns real kind dp dimension pointer rns real kind dp dimension pointer nodes real kind dp dimension 0 2 gO intege
43. routine pstab2n x y px py typ ptyp tabi tab2 real kind dp intent in x y integer intent in px py typ ptyp real kind dp dimension 0 intent out tab1 tab2 Evaluate base functions up to degree px py pz at point x y z subroutine ptab3 x y z px py pz typ ptyp tab real kind dp intent in x y z integer intent in px py pz typ ptyp real kind dp intent out tab 0 Evaluate derivatives of base functions up to degree px py pz at point x y z tabl contains the x derivative tab2 the y derivative and tab3 the z derivative subroutine pstab3 x y z px py pz typ ptyp tabl tab2 tab3 real kind dp intent in x y z integer intent in px py pz typ ptyp real kind dp dimension 0 intent out tab1 tab2 tab3 Interpolates fi xi with n 1 nodes using monomials subroutine ipol1 fi xi n coeff integer intent in n real kind dp intent in fi 0 n xi 0 n real kind dp intent out coeff 0 n Compute all Lagrange base functions with polynomial degrees px on a rectangle or triangle using a tensor product node distribution 24 subroutine ipol2 fi xil xi2 typ px coeff integer intent in px 0 1 typ real kind dp dimension 0 intent in xil xi2 real kind dp dimension 0 mo 0 mo intent in fi real kind dp dimension 0 mo 0 mo intent out coeff Compute the coefficients of the i th Lagrange polynomial with n 1 nodes subroutine lagrangei xi n i lag
44. s integer dimension pointer rnodes integer ng bmode 0 1 character len 32 name type elem2 dimension pointer elem Array of mesh elements allowing hanging nodes integer nelem end type mesh2 rd2 i rp 1 i rx i rdl i rp 0 i rdl i parallelogram rt i 4 triangle rt i 3 4 3 Three dimensional spaces type elem3 integer rnodes 0 8 rnodes 0 8 hexaeder 4 tetraeder rnodes 1 8 node numbers integer rc 0 5 neighbours integer rct 0 5 sides of neighbours integer father integer sons 1 8 integer ref end type elem3 type mesh3 real kind dp dimension pointer nodes integer nnodes integer dimension pointer rnodes integer ng bmode 0 1 integer dimension pointer rc rct character len 32 name type elem3 dimension pointer elem Array of mesh elements allowing hanging nodes integer nelem end type mesh3 4 4 Mesh databases 5 Internal data structures for Splines There are three data types created for the description of splines defined on a manifold with 1 2 or 3 dimensions Note that the dimension of the manifold can be different from the dimension of the space Depending on this dimension the definitions are contained in the modules splinel spline2 or spline3 located in the file splines1 f90 splines2 f90 or splines3 f90 mesh and spline definitions are located in the sa
45. s and their derivatives are given in the modules funct located in the files funct f90 where denotes the short cut of the configuration Here we have to distinguish mainly between Dirichlet and Neumann boundary data and volume data with the prefixes d n and v respectively The operators are denoted by the prefixes lap lame helm and maxw The data function is denoted by the postfix f and the corresponding exact solution is denoted by fe An additional derivative is indicated by a s E g the derivative of Dirichlet data for the Laplacian is denoted by dslapf d s lap f and the exact solution of the Dirichlet problem by dlapfe Additionally there is the module ltyp defined containing the integer variable typ 0 1 steering the choice of the function which has to be evaluated The subroutines implementing the boundary data functions have the following structure subroutine pre op f x nx fx use ltyp real kind dp x 0 nx 0 fx 0 x select case 1typ 0 case 0 case 1 end select end subroutine pre op f Volume data functions are given by subroutine v op f x fx use ltyp real kind dp x 0 x fx 0 x select case ltyp 0 case 0 case 1 end select end subroutine v op f The actual components of the arguments x and nx depends on the configuration nx describes the direction of a derivative usually the normal direction
46. sing the prefix lap lame or helm We usually don t state the whole boundary but we compute all integrals for two boundary elements from which all other terms can be inferred There are two ways of describing a boundary element T First we can state the start point a and the end point e of T or we can state the mid point m e a 2 and the difference vector from end point and mid point d e a 2 Most operators demand the definition of some constants In the case of the Lam equation the Lam coefficients are defined in the program by using subroutine initlame lambda mu real kind dp intent in lambda mu 27 before any other subroutine is called All internal constants which are needed by liblame2 f90 are defined by this subroutine As long as the coefficients don t change it is sufficient to call the subroutine only once In the case of the Helmholtz equation the wavenumber kw and the number of terms of the series expansion kn are defined by subroutine inithelm kw kn real kind dp intent in kw integer intent in kn The test and ansatz functions are given by k fer k 2 q2 cel 8 otherwise The Galerkin elements are usually given by Fri ote mA dde To Fi II 1 1 Idol al f J tht F do ty mo dy ty m1 dty dt 241 The values of the Galerkin elements can be real numbers complex numbers or matrices They are always sto
47. stribution of the nested subspaces is given by integer mtop am 0 mtop xfm 0 mtop xfn 0 mtop The projections restriction and prolongation from level to level are given by the following arrays integer prom parameter prom 9 integer pro 0 prom 0 hpro 0 prom 0 real prow 1 prom 0 hprow 1 prom 0 integer difsi 0 5 0 real difsw 1 5 0 Using these definitions we implement the Multilevel additive Schwarz method with exact solution at the lowest level subroutine bmasv2 ktyp mtop mdc am a xfm xfn x f pro prow ix difsi difsw rp ri ngm ngn hpro hprow haard ipvt c computes one mas step at level mtop c exact solution at lowest level 42 integer ktyp mdc real a 0O 0 x 0 real haard 0 integer ix 0 ipvt 0 integer rp 0 ri 0 1 0 c Position und Anzahl der Netze integer ngm 0 mtop ngn 0 mtop the multigrid algorithm subroutine mg mtop mu nul nu2 mds mdc mgr omega am a xfm xfn x f pro prow ix difsi difsw drh rp ri rci ngm ngn ncm ncn dntyp ktyp nt ipvt c fuehrt einen Multigridschritt auf level mtop durch integer mu dntyp ktyp nt mds mdc mgr integer nul 0 mtop nu2 0 mtop real omega 0 mtop real a 0 f 0 x 0 integer ix 0 rci 0 ipvt 0 integer rp O nt 1 0 ri 0 1 0 real drh 0 2 0 c Position und Anzahl der Netze integer ngm 0 mtop ngn 0 mtop c Position und Anzahl der logischen Freiheitsgrade integer ncm 0 mtop ncn 0 mt
48. t in ckom 0 x 0 2 type spline3 intent in sp real kind dp intent in optional a 0 2 0 7 real kind dp intent out fwx 0 fwy 0 fwz 0 5 4 Spline databases There exist three databases for spline definitions which contain all managed spline definitions for 1 d 2 d and 3 d The databases can be queried by the following pointer valued functions qspli qspl2 qspl3 e g type spline1 pointer spl spi gt qspl1 Name of Spline definition if not associated sp1 then write msg Spline definition not found call logmsg end if The coefficients of representations of this spline definitions can be queried by qsp11x qsp12x qsp13x This results to pointers to a data structure which contains a pointer to the spline definition and a pointer to the vector of coefficients Note that there can be more than one spline using a single spline definition type spline1x pointer spix type spline1 pointer spl real kind dp dimension pointer x spix gt qsplix Name of Spline if associated spix then spixfsp is pointer to spline definition spix x is pointer to coefficient vector sp1 gt sp1x sp x gt sp1x x end if 12 5 5 bcdat database Every spline class has its own local database for additional information which will not be stored in its own component of the Fortran data structure It is not necessary to change fundamental data structures if you simply want to in
49. the arc segment whereas the last node defines a corner 2 4 curved triangular element sphere shell The first three nodes define the cor ners of the curved triangle the last node is a node in the curved triangle which allows to determine the center of the sphere 2 6 defines a curved quadrilateral element cylinder shell The first three nodes define the first arc segment the last three nodes the second arg segment 10 n2 Curvi linear element square represented by Lagrange based uniform interpolation polynomial with n nodes bmode 3 3 type nr 0 4 tetrahedral element 0 5 pyramidal element the first four nodes give the base of the pyramid 0 6 prismatic element the first three nodes define the bottom 0 8 hexahedral element the first four nodes define the bottom 1 8 cylinder with one curved surface The first 4 nodes define one end of the cylinder the next 4 nodes the other end 1 10 cylinder this time based on a quadrilateral element 1 12 hollow cylinder shell not fully working yet The first 6 nodes define one end of the cylinder the second 6 nodes the other end From each set of 6 nodes the first 3 nodes describe the outer arc and the other ones the second arc cf bmode 2 2 typ 1 nr 6 2 4 sphere The first node is the center the last three nodes define the surface segment 2 8 spherical segment with two curved surfaces The first three nodes define the corners of the first curved tria
50. tinguish two cases i e two reference elements First we have the rectangle Q 1 1 Here we can use tensor products of antiderivatives of Legendre polynomials as base functions i e ralz y Le e Lily OSk lt pe OSI lt py If we take at look at this in matrix form we obtain Opy t 01 01 11 00 10 0 gt p 0 We can read this as follows poo lower left corner Pio lower right corner Qoi upper left corner P10 upper right corner ko k gt 2 lower edge functions k k gt 2 upper edge functions go 1 gt 2 left edge functions du 1 gt 2 right edge functions Phi inner bubble functions Second we have the triangle T xz y 0 lt yy S 1 z 0 lt z lt 1 Here we have the difficulties that we don t use the interval 0 1 instead of the standard interval 1 1 as before therefore we have to use the transformed polynomials Ly z L 2x 1 x 0 1 and the main difficulty We can not use tensor products anymore to construct continuous splines But we have to ensure that our splines at triangles can be connected easily with splines defined at rectangles We achieve this goal by using the polynomials Ly 1 L x 1 z k gt 2 This leads to the definitions 19 goo z y 1 y lower left corner diollz y z lower right corner boi z y y upper left corner Prolz y Le x yLe x 2 lt k lt p lower edge polz y Lily aLl
51. tions are a linear combination of local basis functions Every local basis functions belongs only to one and only one global basis function Using this fundamental assumptions we have the following objects to deal with e The mesh consisting of mesh elements w e The mappings F Q gt wi from the reference element to the mesh elements e The set of basis function on the reference element Q gt R or R3 e The local basis functions on every element ne z A EA 2 wi R or if pref is vector valued 61 x BPE 2 w gt R with a transformation matrix B RIXA Then every global basis function can be represented in the following way Ppl 2 wi edie a 1 p ri k Here w are weights belonging to every local basis functions and r denotes the global basis functions to which every local basis function belongs to First we start with the elements which are in common for all cases type spline 1 2 3 integer ng ktyp vtypc vtypf dof integer dimension pointer ri gt null ri 0 ng integer dimension pointer rci gt null rci 0 real kind dp dimension pointer rcw gt null rcw 0 ktyp 1 means we deal with real numbers and ktyp 2 means we are dealing with complex valued numbers vtypc denotes the number of coefficients belonging to a basis function vtypf 1 means we deal with a solution with scalar components whereas vtypf 2 3 means we are dealing with a vector valued
52. troduce new information Entries in the database are created by the generic subroutines crbcdat and addbcdat crbcdat simply creates the entry in the database using the fullname and nickname and reserves the amount of memory given by ni integer nr real or nc character addbcdat uses directly the given data i r c subroutine crbcdat sp fullname nickname ni nr nc type spline intent inout sp character len intent in fullname nickname integer intent in optional ni nr nc subroutine addbcdat sp fullname nickname i r c type spline intent inout sp character len intent in fullname nickname integer intent in optional i real kind dp intent in optional r character len intent in optional c Using the name fullname or nickname you can query the bedat database in the usual way Namely by qbci for integer data by qbcr for real data and by qbcc for strings function qbci sp name integer dimension pointer qbci type spline intent in sp character len intent in name function qbcr sp name real kind dp dimension pointer qbcr type spline intent in sp character len intent in name function qbcc sp name type string pointer qbcc type spline intent in sp character len intent in name 6 Problem data structure 6 1 matstr Using matstr we describe the construction of the block system matrix Entries ar
53. which are orthonormal with respect to this moments Let Jf tq ds if q pr e mi o 4 fpo x nqdf ifq Qk 2 k 1 X Qk 1 k 2 7 Job q da if q Qk 1 k 2 k 2 X Qk 2 k 1 k 2 X Qk 2 k 2 k 1 Then the basis functions i have to satisfy the condition m5 bi We can represent 6 with respect to di Le 3k k 1 i 5 aip l 1 This leads to the following linear system 3k k41 m di 5 aim thr ds l 1 If we define bj m 15 then we have aa a bj 50 22 10 2 4 Nedelec elements N D on tetrahedrals 11 Polynomials Most of the subroutines for use and manipulation of polynomials are contained in the mod ule poly located in the file fox poly f90 The module polydat contains some global data structures created to avoid the recomputation of polynomial coefficients Most important is the definition of the polynomial basis functions They are declared by the integer parameter ptyp In the moment the following values for ptyp are allowed typ basis Monomials Legendre polynomials Tschebyscheff polynomials Antiderivatives of Legendre polynomials Orthonormal Legendre polynomials Lagrange polynomials uniform nodes 6 Lagrange polynomials arbitrary nodes 10 Raviart Thomas RT polynomials 13 Nedelec ND polynomials 14 Traces of Nedelec TND polynomials 15 PEERS polynomials oP Nr OD The type of polynomials for a given ptyp can be printed in clear t
54. y ty amp pvx pvy pkx pky pwx pwy amp amp vklmn kklmn wklmn Element type 3 triangle 4 rectangle integer intent in tx ty vector of the origin boundaries normal direction real kind dp dimension 0 2 intent in bx dx1 dx2 nx real kind dp dimension 0 2 intent in by dy1 dy2 ny polynomial degree integer dimension 0 1 intent in px pkx pwx integer dimension 0 1 intent in py pky pwy Galerkin matrices real kind dp dimension 0 fm 0 fm 0 fm 0 fm 0 nt 1 intent inout vklmn kklmn wklmn In the same way we can compute the Galerkin elements for the adjoint double layer potential which are usually computed by taking the transposed matrix of the double layer potential subroutine integ3ks bx dx1 dx2 nx tx by dy1 dy2 ny ty amp amp pkx pky ksklmn Element type 3 triangle 4 rectangle integer intent in tx ty vector of the origin boundaries normal direction real kind dp dimension 0 2 intent in bx dx1 dx2 nx real kind dp dimension 0 2 intent in by dy1 dy2 ny polynomial degree integer intent in pkx 0 1 pky 0 1 Galerkin matrix real kind dp intent inout ksklmn 0 fm 0 fm 0 fm 0 fm 0 nt 1 For the computation of the right hand side we need the integrals over a general function 34 multiplied with the test functions Ta F Z po ki T ds To subroutine genid23 bx dx1 dx2 nx tx px gqn gqk gqw kid fkt nt Element typ 3 t
55. z y ne gt llor Ven a cry TET gt un Alien TETh More ln Es z2 r y DW2F u ETAF DW2JSIGMA u EJSIG RES2CXIMTXI u D N XIMTXI 1 K N V p n ne 122 E Erllz2c ew RES2CPTXIS p D N PTXIS 1 K N V p 1m Mpe A pn t Ef lleen RES2CXIMTXIF p S RES2CPTXISF p S RES2CPTFREEF p S RES2CZETA u D N ZETA C I K p n WD ne Ve hellSllZ2 e LSQ2CDIVPF p DIVPF div pa flirer LSQ2CNABLA u p NABLA p Vunlif2 r LSQ2CVARPHI u D N LSQ2CDIFXI u D N DIFXI W un uo 2pn n 2to I K on to Ils 2 6 PLAP2CETAGR u PLAP2CETAF u RES2CTRANSX u D N 18 10 Local basis functions 10 1 Scalar basis functions The routines for manipulating polynomials are described in their own section Here we deal with the local properties of continuous splines and their representation Mainly we use antiderivatives of Legendre polynomials to ensure the continuity In one dimension this is not a problem at all we just use the definition of the antiderivatives of Legendre polynomials at a reference element e g Lola 1 0 2 A z 1 z 2 Ln a La z In 2 2n 1 YE Ln 1 y dy where L x is a Legendre polynomial defined by Lola 1 L z z n D Enyi a 2n 1 z L z nIn z The antiderivatives of Legendre polynomials have the important and simplifying property L n 1 0 n gt 2 In two dimensions we have to dis
Download Pdf Manuals
Related Search