The Mathematica Journal
Departments
Feature Articles
Columns
New Products
New Publications
Classifieds
Calendar
News Bulletins
Mailbox
Letters
Write Us
About the Journal
Staff and Contributors
Submissions
Subscriptions
Advertising
Back Issues
Home
Download this Issue

Going Global

Element-level operations establish the local equations governing the behavior of the element in the subdomain [Graphics:../Images/index_gr_81.gif] of the element. However, the behavior of the system, which is the assembly of all the elements, is defined on the domain [Graphics:../Images/index_gr_82.gif], where is the assembly operator.

In practical terms this amounts to an assembly of the local element integrals into a global system of equations. The assembly operation is best described with an example. Let's create a mesh with a set of LaplacianElements as described in the previous section.

[Graphics:../Images/index_gr_83.gif]

Three rows and four columns of nodes describe this mesh. These nodes are created as a regular grid with the following statements.

[Graphics:../Images/index_gr_84.gif]

Next, two rows and three columns of elements on the regular grid are described as:

        els = Flatten[Table[ 
            new[LaplacianElement,
                {nods[[i +   (j-1)*cols]], 
                 nods[[i + 1+(j-1)*cols]],
                 nods[[i + 1+    j*cols]], 
                 nods[[i +       j*cols]]}]
            ,{j,1,rows-1}, {i,1,cols-1} ]] ;

Let's create a mesh describing this geometry and the topology of elements. We have created 12 nodes and 6 elements in two lists called nods and els respectively. The mesh is an instance of the class called Fem.

        mesh = new[Fem, nods, els];
        Show[ display[mesh], ViewPoint->{-7.091, -2.682, 4.449}];

[Graphics:../Images/index_gr_85.gif]

Using the element-level matrices, the global system of equations can be described for the Laplacian equation on this grid. This involves computing the element integrals and assembling them into the global system of equations. The assembly operation is performed by the method called assembleblock in the Fem class. This method requires a matrix-valued element function such as Laplacian[e_] which was defined above. It will then loop over the elements in the mesh and assemble the matrices returned by the function passed in as the argument into the global matrix, K. In this example the global matrix K is a [Graphics:../Images/index_gr_86.gif] matrix.

        K = assembleblock[mesh,Laplacian,1];
        MatrixForm[K]

The Laplacian field on this grid can be calculated for a given set of boundary conditions. If the boundary conditions are defined as [Graphics:../Images/index_gr_88.gif] at [Graphics:../Images/index_gr_89.gif] and [Graphics:../Images/index_gr_90.gif], the exact solution of this problem exists and is defined by [Graphics:../Images/index_gr_91.gif].

The boundary conditions describe the values of the Laplacian field at the nodes 1, 5, 9 and 4, 8, 12. The values of the Laplacian field for the remainder of the nodes must be computed. The global system of equations is given as:

[Graphics:../Images/index_gr_92.gif]

The condensed form of the equations, after the prescribed values multiply their columns, are moved to the right-hand side of the equations, and written as:

[Graphics:../Images/index_gr_93.gif]

The resulting equations can be solved with:

    nodalvalues  = {{0},{0},{0},{3},{0},{0},{0},{3},{0},{0},{0},{3}};
    unknownnodes = {     2,  3,          6,  7,          10, 11};
    Kcondensed = K[[unknownnodes ,unknownnodes ]];
    load       = (K  . nodalvalues )[[unknownnodes ,1]];
    LinearSolve[Kcondensed,-load]
[Graphics:../Images/index_gr_94.gif]

These are the exact answers for the nodes {2, 3, 6, 7, 10, 11} respectively. However, the Fem package provides another convenience for imposing boundary conditions and solving the resulting linear system. I can redefine the same problem as:

        fixednodes       = {{ 1,4,5,8,9,12}};
        fixedvalues      = {{ 0,3,0,3,0,3 }};
        b                = Table[ 0, {i,1,12}, {j,1,1} ];
        linearSolve[mesh,K,b,fixednodes,fixedvalues ]
[Graphics:../Images/index_gr_95.gif]


Converted by Mathematica    

[Article Index] [Prev Page][Next Page]