A simple FEM Code

Version 0.1, April 2000
Andreas Stahel, HTA Biel, Switzerland
Andreas.Stahel@hta-bi.bfh.ch

Description

This package allows to solve second order elliptic differential equations in two variables:

     div(a*grad u) - b*u = f         in domain
                                      u =  gD     Dirichlet boundary conditions on first part of boundary
                         a*du/dn =  gN      Neumann condition on the other part of the boundary
                       
If the functions a, b f, gD and gN are given, then a numerical approximation is computed, using the method of finite elements.

The mesh on the domain has to be generated by the program EasyMesh.  Within EasyMesh the boundary points have to be marked, such that the code knows whether to apply Dirichlet or Neumann boundary conditions. The codes are
      type 1  = Dirichlet boundary
      type 2  = Neumann boundary
The files EasyMeshIntro.A4.ps and EasyMeshIntro.letter.ps may help to get started. But as usual the standard documentation is best, see at http://www-dinma.univ.trieste.it/~nirftc/research/easymesh/ 
If you plan to use this package you will have to install EasyMesh. The enclosed examples can be run without installing the mesh generator.

The command given in this package should
-  read the information on the mesh.
-  create the finite element model, set up the system of linear equations and solve it.
-  allow to visualize the solutions.
Consider the enclosed example.

The notebook FEM2DExamples.nb contains a few examples.
The notebook boundarygenerator.nb ( written by Pierre Andre Chevalier) does allow to generate the input file for EasyMesh within a Mathematica notebook.
The subdirectory ./EasyMesh contains some input files for the examples.


This code was written as support material for a class on Calculus of Variations and Finite Elements. It should only be used for educational purposes.  It was developed and tested on a Linux system and thus will 
behave the same on all Unix systems. I will try to help to get it up and running. On MS basesd systems you are on your own.
- If you like it, let me know
- If you have suggestions on how to impove the package, let me know. 
   There is some obvious space for improvement.
- If you think it is a piece of .... , delete it and dont tell anybody 

Enjoy this day

Andreas

A First Example

Reading the Package and display the built in documentation

<< BVP2.m

?ReadMesh

?MeshGraphics

?FEMSolve

?FEMGradient

?FEMSurface

?FEMSurfaceNoMesh

?LevelGraphics2D

?LevelGraphics3D

Visualization of Mesh for Control

{nodes, elements, segments} = ReadMesh["./example2"];

gmesh = MeshGraphics[nodes, elements];

\(\(Show[gmesh\ , 
AspectRatio -> Automatic, \ 
Axes -> Automatic]; 
    \)\)

Computations

gD[x_, y_] = 
      Which[3.5 ≤ x ≤ 4.5 && 
          0.5 ≤ y ≤ 1.5, 100, True, 0]; a = Function[{x, y}, 1]; b = Function[{x, y}, 0]; f = Function[{x, y}, 0]; gDirichlet = Function[{x, y}, gD[x, y]]; gNeumann = Function[{x, y}, 1]; points = 
      FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, segments];

A 3D Plot

g = FEMSurface[points, elements]; \(\(solution = 
      Show[g, 
		Axes -> True, 
        AxesLabel -> {"x", "y", "u"}, AspectRatio -> 1, 
        PlotRange -> All, ViewPoint -> {1.582, \(-2.9\), 0.647}]; 
    \)\)

Level curves

\(levels = Table[j, {j, 0, 1, 0.05}]; 

    glevel = LevelGraphics2D[points, elements, levels]; 
Show[glevel]; \)

\(levels = Table[j, {j, 0, 1, 0.05}]; 

    glevel3D = LevelGraphics3D[points, elements, levels]; 
Show[glevel3D]; 
    \)

\(g = FEMSurfaceNoMesh[points, elements]; 

    solution = 
      Show[g, 
		Axes -> True, 
        AxesLabel -> {"x", "y", "u"}, AspectRatio -> 1, 
        PlotRange -> All, ViewPoint -> {1.582, \(-2.9\), 0.647}]; 
    Show[solution, glevel3D]; \)

Internal Documentation

Global variables and functions

The code in the package is using the global variables nodes, elements, segments and Node2Degree.
Modifying those may cause problems and lead to wrong results.

Internal variables and functions

{dof , Node2Degree} = FindDOF[nodes]
   dof is the number of degrees of freedom
   Node2Degree contains the mapping of the nodes to the DOF, an entry 0 corresponds to a fixed value
   
 {elMat, elVec} =  ElementContribution[corners, aCoeff, bCoeff, fCoeff]
    find the  element stiffness matrix and vector for an element with given corners. 
    aCoeff, bCoeff, fCoeff are the functions describing the PDE 
   
segVec=ElementContributionEdge[corners, gCoeff]
    find the contribution of one segment, Nemann condition with function gCoeff   
   
insertElement[element, aCoeff, bCoeff, fCoeff, gCoeff]
    insert the conribution of the element with the given nodes into Aglobal and bglobal
    
insertElementEdge[sdge, gCoeff]
    insert the conribution of the element with the given nodes into  bglobal
   
 Aglobal  the global stiffness matrix
 bGlobal  the global vector