~chaffra/fiat/main

« back to all changes in this revision

Viewing changes to FIAT/Lagrange.py

  • Committer: kirby
  • Date: 2005-05-31 15:03:57 UTC
  • Revision ID: devnull@localhost-20050531150357-n7ka5987yo6ylpaz
[project @ 2005-05-31 15:03:57 by kirby]
Initial revision

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
# Written by Robert C. Kirby
 
2
# Copyright 2005 by The University of Chicago
 
3
# Distributed under the LGPL license
 
4
# This work is partially supported by the US Department of Energy
 
5
# under award number DE-FG02-04ER25650
 
6
 
 
7
# Last modified 16 May 2005
 
8
 
 
9
 
 
10
"""This module defines the standard Lagrange finite element over
 
11
any shape for which an appropriate expansion basis is defined.  The
 
12
module simply provides the class Lagrange, which takes a shape and a
 
13
degree >= 1.
 
14
"""
 
15
 
 
16
import dualbasis, functional, points, shapes, polynomial, functionalset
 
17
 
 
18
class LagrangeDual( dualbasis.DualBasis ):
 
19
    """Dual basis for the plain vanilla Lagrange finite element."""
 
20
    def __init__( self , shape , n , U ):
 
21
        # Nodes are point evaluation on the lattice, which here is
 
22
        # ordered by topological dimension and entity within that
 
23
        # dimension.  For triangles we have vertex points plus
 
24
        # n - 1 points per edge.  We then have max(0,(n-2)(n-1)/2)
 
25
        # points on the interior of the triangle, or that many points
 
26
        # per face of tetrahedra plus max(0,(n-2)(n-1)n/6 points
 
27
        # in the interior
 
28
        pts = [ [ points.make_points(shape,d,e,n) \
 
29
                      for e in shapes.entity_range(shape,d) ] \
 
30
                    for d in shapes.dimension_range(shape) ]
 
31
 
 
32
        pts_flat = reduce( lambda a,b:a+b , \
 
33
                           reduce( lambda a,b:a+b , pts ) )
 
34
 
 
35
        # get node location for each node i by self.pts[i]
 
36
        self.pts = pts_flat
 
37
 
 
38
        ls = functional.make_point_evaluations( U , pts_flat )
 
39
 
 
40
        # entity_ids is a dictionary whose keys are the topological
 
41
        # dimensions (0,1,2,3) and values are dictionaries mapping
 
42
        # the ids of entities of that dimension to the list of
 
43
        # ids of nodes associated with that entity.
 
44
        # see FIAT.base.dualbasis for description of entity_ids
 
45
        entity_ids = {}
 
46
        id_cur = 0
 
47
        for d in shapes.dimension_range( shape ):
 
48
            entity_ids[d] = {}
 
49
            for v in shapes.entity_range( shape, d ):
 
50
                num_nods = len( pts[d][v] )
 
51
                entity_ids[d][v] = range(id_cur,id_cur + num_nods)
 
52
                id_cur += num_nods
 
53
 
 
54
        dualbasis.DualBasis.__init__( self , \
 
55
                                      functionalset.FunctionalSet( U , ls ) , \
 
56
                                      entity_ids )
 
57
                                      
 
58
        return
 
59
 
 
60
class Lagrange( polynomial.FiniteElement ):
 
61
    def __init__( self , shape , n ):
 
62
        if n < 1:
 
63
            raise RuntimeError, \
 
64
                  "Lagrange elements are only defined for n >= 1"
 
65
        U = polynomial.OrthogonalPolynomialSet( shape , n )
 
66
        Udual = LagrangeDual( shape , n , U )
 
67
        polynomial.FiniteElement.__init__( self , \
 
68
                                           Udual , U )
 
69
 
 
70
class VectorLagrangeDual( dualbasis.DualBasis ):
 
71
    def __init__( self , shape , n , U ):
 
72
        space_dim = shapes.dimension( shape )
 
73
        nc = U.tensor_shape()[0]
 
74
        pts = [ [ points.make_points(shape,d,e,n) \
 
75
                  for e in shapes.entity_range(shape,d) ] \
 
76
                for d in shapes.dimension_range(shape) ]
 
77
 
 
78
        pts_flat = reduce( lambda a,b:a+b , \
 
79
                           reduce( lambda a,b:a+b , pts ) )
 
80
 
 
81
        self.pts = pts_flat
 
82
        
 
83
        ls = [ functional.ComponentPointEvaluation( U , c , pt ) \
 
84
               for c in range(nc) \
 
85
               for pt in pts_flat ]
 
86
        
 
87
 
 
88
        # see FIAT.base.dualbasis for description of entity_ids
 
89
        entity_ids = {}
 
90
        id_cur = 0
 
91
        for d in shapes.dimension_range( shape ):
 
92
            entity_ids[d] = {}
 
93
            for v in shapes.entity_range( shape, d ):
 
94
                num_nods = len( pts[d][v] )
 
95
                entity_ids[d][v] = range(id_cur,id_cur + num_nods)
 
96
                id_cur += num_nods
 
97
        
 
98
        fset = functionalset.FunctionalSet( U , ls )
 
99
 
 
100
        dualbasis.DualBasis.__init__( self , fset , entity_ids , nc ) 
 
101
 
 
102
 
 
103
class VectorLagrange( polynomial.FiniteElement ):
 
104
    def __init__( self , shape , n , nc = None ):
 
105
        if n < 1:
 
106
            raise RuntimeError, \
 
107
                  "Lagrange elements are only defined for n >= 1"
 
108
        U = polynomial.OrthogonalPolynomialArraySet( shape , n , nc )
 
109
        Udual = VectorLagrangeDual( shape , n , U )
 
110
        polynomial.FiniteElement.__init__( self , Udual , U )