← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/fem/dolfin/StiffnessMatrix2D.h

  • Committer: Anders Logg
  • Date: 2007-01-10 09:04:44 UTC
  • mfrom: (1689.1.221 trunk)
  • Revision ID: logg@simula.no-20070110090444-ecyux3n1qnei4i8h
Remove old head

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/fem/dolfin/StiffnessMatrix2D.h
1
 
// Automatically generated by FFC, the FEniCS Form Compiler, version 0.2.5.
2
 
// For further information, go to http://www/fenics.org/ffc/.
3
 
// Licensed under the GNU GPL Version 2.
4
 
 
5
 
#ifndef __STIFFNESSMATRIX2D_H
6
 
#define __STIFFNESSMATRIX2D_H
7
 
 
8
 
#include <dolfin/Mesh.h>
9
 
#include <dolfin/Cell.h>
10
 
#include <dolfin/Point.h>
11
 
#include <dolfin/Vector.h>
12
 
#include <dolfin/AffineMap.h>
13
 
#include <dolfin/FiniteElement.h>
14
 
#include <dolfin/FiniteElementSpec.h>
15
 
#include <dolfin/LinearForm.h>
16
 
#include <dolfin/BilinearForm.h>
17
 
 
18
 
namespace dolfin { namespace StiffnessMatrix2D {
19
 
 
20
 
/// This class contains the form to be evaluated, including
21
 
/// contributions from the interior and boundary of the domain.
22
 
 
23
 
class BilinearForm : public dolfin::BilinearForm
24
 
{
25
 
public:
26
 
  
27
 
  class TestElement : public dolfin::FiniteElement
28
 
  {
29
 
  public:
30
 
  
31
 
    TestElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)
32
 
    {
33
 
      // Element is scalar, don't need to initialize tensordims
34
 
  
35
 
      // Element is simple, don't need to initialize subelements
36
 
    }
37
 
  
38
 
    ~TestElement()
39
 
    {
40
 
      if ( tensordims ) delete [] tensordims;
41
 
      if ( subelements )
42
 
      {
43
 
        for (unsigned int i = 0; i < elementdim(); i++)
44
 
          delete subelements[i];
45
 
        delete [] subelements;
46
 
      }
47
 
    }
48
 
  
49
 
    inline unsigned int spacedim() const
50
 
    {
51
 
      return 3;
52
 
    }
53
 
  
54
 
    inline unsigned int shapedim() const
55
 
    {
56
 
      return 2;
57
 
    }
58
 
  
59
 
    inline unsigned int tensordim(unsigned int i) const
60
 
    {
61
 
      dolfin_error("Element is scalar.");
62
 
      return 0;
63
 
    }
64
 
  
65
 
    inline unsigned int elementdim() const
66
 
    {
67
 
      return 1;
68
 
    }
69
 
  
70
 
    inline unsigned int rank() const
71
 
    {
72
 
      return 0;
73
 
    }
74
 
  
75
 
    void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
76
 
    {
77
 
      nodes[0] = cell.vertexID(0);
78
 
      nodes[1] = cell.vertexID(1);
79
 
      nodes[2] = cell.vertexID(2);
80
 
    }
81
 
  
82
 
    void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
83
 
    {
84
 
      points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);
85
 
      points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);
86
 
      points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);
87
 
      components[0] = 0;
88
 
      components[1] = 0;
89
 
      components[2] = 0;
90
 
    }
91
 
  
92
 
    void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const
93
 
    {
94
 
      // FIXME: Temporary fix for Lagrange elements
95
 
      values[0] = x[vertex];
96
 
    }
97
 
  
98
 
    const FiniteElement& operator[] (unsigned int i) const
99
 
    {
100
 
      return *this;
101
 
    }
102
 
  
103
 
    FiniteElement& operator[] (unsigned int i)
104
 
    {
105
 
      return *this;
106
 
    }
107
 
  
108
 
    FiniteElementSpec spec() const
109
 
    {
110
 
      FiniteElementSpec s("Lagrange", "triangle", 1);
111
 
      return s;
112
 
    }
113
 
    
114
 
  private:
115
 
  
116
 
    unsigned int* tensordims;
117
 
    FiniteElement** subelements;
118
 
  
119
 
  };
120
 
    
121
 
  class TrialElement : public dolfin::FiniteElement
122
 
  {
123
 
  public:
124
 
  
125
 
    TrialElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)
126
 
    {
127
 
      // Element is scalar, don't need to initialize tensordims
128
 
  
129
 
      // Element is simple, don't need to initialize subelements
130
 
    }
131
 
  
132
 
    ~TrialElement()
133
 
    {
134
 
      if ( tensordims ) delete [] tensordims;
135
 
      if ( subelements )
136
 
      {
137
 
        for (unsigned int i = 0; i < elementdim(); i++)
138
 
          delete subelements[i];
139
 
        delete [] subelements;
140
 
      }
141
 
    }
142
 
  
143
 
    inline unsigned int spacedim() const
144
 
    {
145
 
      return 3;
146
 
    }
147
 
  
148
 
    inline unsigned int shapedim() const
149
 
    {
150
 
      return 2;
151
 
    }
152
 
  
153
 
    inline unsigned int tensordim(unsigned int i) const
154
 
    {
155
 
      dolfin_error("Element is scalar.");
156
 
      return 0;
157
 
    }
158
 
  
159
 
    inline unsigned int elementdim() const
160
 
    {
161
 
      return 1;
162
 
    }
163
 
  
164
 
    inline unsigned int rank() const
165
 
    {
166
 
      return 0;
167
 
    }
168
 
  
169
 
    void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
170
 
    {
171
 
      nodes[0] = cell.vertexID(0);
172
 
      nodes[1] = cell.vertexID(1);
173
 
      nodes[2] = cell.vertexID(2);
174
 
    }
175
 
  
176
 
    void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
177
 
    {
178
 
      points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);
179
 
      points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);
180
 
      points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);
181
 
      components[0] = 0;
182
 
      components[1] = 0;
183
 
      components[2] = 0;
184
 
    }
185
 
  
186
 
    void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const
187
 
    {
188
 
      // FIXME: Temporary fix for Lagrange elements
189
 
      values[0] = x[vertex];
190
 
    }
191
 
  
192
 
    const FiniteElement& operator[] (unsigned int i) const
193
 
    {
194
 
      return *this;
195
 
    }
196
 
  
197
 
    FiniteElement& operator[] (unsigned int i)
198
 
    {
199
 
      return *this;
200
 
    }
201
 
  
202
 
    FiniteElementSpec spec() const
203
 
    {
204
 
      FiniteElementSpec s("Lagrange", "triangle", 1);
205
 
      return s;
206
 
    }
207
 
    
208
 
  private:
209
 
  
210
 
    unsigned int* tensordims;
211
 
    FiniteElement** subelements;
212
 
  
213
 
  };
214
 
  
215
 
  BilinearForm(const real& c0) : dolfin::BilinearForm(0), c0(c0)
216
 
  {
217
 
    // Create finite element for test space
218
 
    _test = new TestElement();
219
 
 
220
 
    // Create finite element for trial space
221
 
    _trial = new TrialElement();
222
 
  }
223
 
 
224
 
  void eval(real block[], const AffineMap& map) const
225
 
  {
226
 
    // Compute geometry tensors
227
 
    const real G0_0_0 = map.det*c0*(map.g00*map.g00 + map.g01*map.g01);
228
 
    const real G0_0_1 = map.det*c0*(map.g00*map.g10 + map.g01*map.g11);
229
 
    const real G0_1_0 = map.det*c0*(map.g10*map.g00 + map.g11*map.g01);
230
 
    const real G0_1_1 = map.det*c0*(map.g10*map.g10 + map.g11*map.g11);
231
 
 
232
 
    // Compute element tensor
233
 
    block[0] = 4.999999999999998e-01*G0_0_0 + 4.999999999999997e-01*G0_0_1 + 4.999999999999997e-01*G0_1_0 + 4.999999999999996e-01*G0_1_1;
234
 
    block[1] = -4.999999999999998e-01*G0_0_0 - 4.999999999999997e-01*G0_1_0;
235
 
    block[2] = -4.999999999999997e-01*G0_0_1 - 4.999999999999996e-01*G0_1_1;
236
 
    block[3] = -4.999999999999998e-01*G0_0_0 - 4.999999999999997e-01*G0_0_1;
237
 
    block[4] = 4.999999999999998e-01*G0_0_0;
238
 
    block[5] = 4.999999999999997e-01*G0_0_1;
239
 
    block[6] = -4.999999999999997e-01*G0_1_0 - 4.999999999999996e-01*G0_1_1;
240
 
    block[7] = 4.999999999999997e-01*G0_1_0;
241
 
    block[8] = 4.999999999999996e-01*G0_1_1;
242
 
  }
243
 
 
244
 
private:
245
 
 
246
 
  const real& c0;
247
 
 
248
 
};
249
 
 
250
 
} }
251
 
 
252
 
#endif