********************************************************************* ******** hplab2D (release 1) ******** ********************************************************************* March 8, 2010 Tomas Vejchodsky www.math.cas.cz/vejchod/hplab.html An hp-FEM solver for Matlab. It solves 2D linear elliptic problems in the form -div( a grad u) + c u = f in Omega, u = 0 on the boundary of Omega by the hp-version of the finite element method (hp-FEM). This version of the FEM is characterized by approximation of the solution u by a piecewise polynomial function with various polynomial degrees on various elements. Proper distribution of polynomial degrees in the mesh enables to achieve exponential rate of convergence even for problems with singularities. The input are the coefficients a and c, the right-hand side f and the hp-FEM mesh in the form of three lists: a list of vertices, a list of triangles, and a list of polynomial degrees. The hp-FEM mesh defines implicitly the polygonal domain Omega. Quick Start =========== - Unpack the archive hplab2Drel1.tar.gz with source files. Keep everything in the subdirectory hplab2D. - Run Matlab. - In Matlab, change current directory to hplab2D. - Run the test computation. >> test1 User's Guide ============ The main function is "hpFEMsolver.m". It has three usual (visible) input parameters: v, t, polydeg and four hidden inputs: m-files: "a_coef.m", "c_coef.m", "f_coef.m", and a global variable "hpFEM_intlambda". The parameters v, t, polydeg describe the hp-FEM mesh. v ... A matrix Nvert-by-2 containing coordinates of the vertices in the triangulation. The value Nvert stands for the number of vertices. t ... A matrix Ntri-by-3 defining the triangles in the mesh. Each row contains three indices to v. polydeg ... A vector Ntri-by-1 with polynomial degrees of the elements. It is so-called interior polynomial degree as opposed to the edge polynomial degree. This can be a single number. Then all element have the same polynomial degree. The m-files "a_coef.m", "c_coef.m", "f_coef.m" are in fact real functions of two variables and they correspond to the coefficients a, c, and f. They must be functions of two input values x and y providing a single real number as an output. Internally, however, the values in barycentres of elements only are evaluated and, hence, the coefficients a, c, and f, are effectively treated as piecewise constant. Finally, the global variable "hpFEM_intlambda" can be ignored on the input and the code will automatically compute the correct values (certain integrals) in the run time. This, however, may lead to a slowdown. It is recommended to initialize the global variable "hpFEM_intlambda" at the beginning of computation by calling "loadprecompint.m" (loading the values from a binary file) or "precompint.m" (computing the values). The output from "hpFEMsolver.m" is a piecewise polynomial function uhp stored in an internal structure (see "hpFEMsolver.m" for more detail). This structure can be used to obtain a value of uhp at any point from the domain Omega by calling "evaluhploc.m" or to plot a graph of uhp by calling "plotuhp.m". For further information see the description directly in the particular m-files. The test example introduced in Quick Start describes well the ability of the code. All data concerning the particular example are stored in a specifice subdirectory - "ex2_circle" in this case. The script "test1.m" just calls "ex2solve.m" from the subdirectory "ex2_circle". The mesh is loaded from the data file "ex2mesh1A.dat" which was created with the help of the Matlab PDE toolbox. The PDE toolbox uses arrays 'p', 'e', and 't' for the description of the triangulation. The code hplab2D needs vertex array v and triangle array t only. They can be obtained from the PDE toolbox arrays as: >> v = p(1:2,:)'; >> t = t(1:3,:)'; Enjoy computation with hplab2D!