← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/modules/poisson/dolfin/Poisson3D.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/modules/poisson/dolfin/Poisson3D.h
1
 
// Automatically generated by FFC, the FEniCS Form Compiler, version 0.2.5.
 
1
// Automatically generated by FFC, the FEniCS Form Compiler, version 0.3.5.
2
2
// For further information, go to http://www/fenics.org/ffc/.
3
3
// Licensed under the GNU GPL Version 2.
4
4
 
8
8
#include <dolfin/Mesh.h>
9
9
#include <dolfin/Cell.h>
10
10
#include <dolfin/Point.h>
11
 
#include <dolfin/Vector.h>
12
11
#include <dolfin/AffineMap.h>
13
12
#include <dolfin/FiniteElement.h>
14
13
#include <dolfin/FiniteElementSpec.h>
 
14
#include <dolfin/BilinearForm.h>
15
15
#include <dolfin/LinearForm.h>
16
 
#include <dolfin/BilinearForm.h>
 
16
#include <dolfin/Functional.h>
 
17
#include <dolfin/FEM.h>
17
18
 
18
19
namespace dolfin { namespace Poisson3D {
19
20
 
23
24
class BilinearForm : public dolfin::BilinearForm
24
25
{
25
26
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 4;
52
 
    }
53
 
  
54
 
    inline unsigned int shapedim() const
55
 
    {
56
 
      return 3;
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
 
      nodes[3] = cell.vertexID(3);
81
 
    }
82
 
  
83
 
    void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
84
 
    {
85
 
      points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
86
 
      points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
87
 
      points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
88
 
      points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
89
 
      components[0] = 0;
90
 
      components[1] = 0;
91
 
      components[2] = 0;
92
 
      components[3] = 0;
93
 
    }
94
 
  
95
 
    void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const
96
 
    {
97
 
      // FIXME: Temporary fix for Lagrange elements
98
 
      values[0] = x[vertex];
99
 
    }
100
 
  
101
 
    const FiniteElement& operator[] (unsigned int i) const
102
 
    {
103
 
      return *this;
104
 
    }
105
 
  
106
 
    FiniteElement& operator[] (unsigned int i)
107
 
    {
108
 
      return *this;
109
 
    }
110
 
  
111
 
    FiniteElementSpec spec() const
112
 
    {
113
 
      FiniteElementSpec s("Lagrange", "tetrahedron", 1);
114
 
      return s;
115
 
    }
116
 
    
117
 
  private:
118
 
  
119
 
    unsigned int* tensordims;
120
 
    FiniteElement** subelements;
121
 
  
122
 
  };
123
 
    
124
 
  class TrialElement : public dolfin::FiniteElement
125
 
  {
126
 
  public:
127
 
  
128
 
    TrialElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)
129
 
    {
130
 
      // Element is scalar, don't need to initialize tensordims
131
 
  
132
 
      // Element is simple, don't need to initialize subelements
133
 
    }
134
 
  
135
 
    ~TrialElement()
136
 
    {
137
 
      if ( tensordims ) delete [] tensordims;
138
 
      if ( subelements )
139
 
      {
140
 
        for (unsigned int i = 0; i < elementdim(); i++)
141
 
          delete subelements[i];
142
 
        delete [] subelements;
143
 
      }
144
 
    }
145
 
  
146
 
    inline unsigned int spacedim() const
147
 
    {
148
 
      return 4;
149
 
    }
150
 
  
151
 
    inline unsigned int shapedim() const
152
 
    {
153
 
      return 3;
154
 
    }
155
 
  
156
 
    inline unsigned int tensordim(unsigned int i) const
157
 
    {
158
 
      dolfin_error("Element is scalar.");
159
 
      return 0;
160
 
    }
161
 
  
162
 
    inline unsigned int elementdim() const
163
 
    {
164
 
      return 1;
165
 
    }
166
 
  
167
 
    inline unsigned int rank() const
168
 
    {
169
 
      return 0;
170
 
    }
171
 
  
172
 
    void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
173
 
    {
174
 
      nodes[0] = cell.vertexID(0);
175
 
      nodes[1] = cell.vertexID(1);
176
 
      nodes[2] = cell.vertexID(2);
177
 
      nodes[3] = cell.vertexID(3);
178
 
    }
179
 
  
180
 
    void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
181
 
    {
182
 
      points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
183
 
      points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
184
 
      points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
185
 
      points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
186
 
      components[0] = 0;
187
 
      components[1] = 0;
188
 
      components[2] = 0;
189
 
      components[3] = 0;
190
 
    }
191
 
  
192
 
    void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const
193
 
    {
194
 
      // FIXME: Temporary fix for Lagrange elements
195
 
      values[0] = x[vertex];
196
 
    }
197
 
  
198
 
    const FiniteElement& operator[] (unsigned int i) const
199
 
    {
200
 
      return *this;
201
 
    }
202
 
  
203
 
    FiniteElement& operator[] (unsigned int i)
204
 
    {
205
 
      return *this;
206
 
    }
207
 
  
208
 
    FiniteElementSpec spec() const
209
 
    {
210
 
      FiniteElementSpec s("Lagrange", "tetrahedron", 1);
211
 
      return s;
212
 
    }
213
 
    
214
 
  private:
215
 
  
216
 
    unsigned int* tensordims;
217
 
    FiniteElement** subelements;
218
 
  
219
 
  };
220
 
  
221
 
  BilinearForm() : dolfin::BilinearForm(0)
222
 
  {
223
 
    // Create finite element for test space
224
 
    _test = new TestElement();
225
 
 
226
 
    // Create finite element for trial space
227
 
    _trial = new TrialElement();
228
 
  }
229
 
 
230
 
  void eval(real block[], const AffineMap& map) const
231
 
  {
232
 
    // Compute geometry tensors
233
 
    const real G0_0_0 = map.det*map.g00*map.g00 + map.det*map.g01*map.g01 + map.det*map.g02*map.g02;
234
 
    const real G0_0_1 = map.det*map.g00*map.g10 + map.det*map.g01*map.g11 + map.det*map.g02*map.g12;
235
 
    const real G0_0_2 = map.det*map.g00*map.g20 + map.det*map.g01*map.g21 + map.det*map.g02*map.g22;
236
 
    const real G0_1_0 = map.det*map.g10*map.g00 + map.det*map.g11*map.g01 + map.det*map.g12*map.g02;
237
 
    const real G0_1_1 = map.det*map.g10*map.g10 + map.det*map.g11*map.g11 + map.det*map.g12*map.g12;
238
 
    const real G0_1_2 = map.det*map.g10*map.g20 + map.det*map.g11*map.g21 + map.det*map.g12*map.g22;
239
 
    const real G0_2_0 = map.det*map.g20*map.g00 + map.det*map.g21*map.g01 + map.det*map.g22*map.g02;
240
 
    const real G0_2_1 = map.det*map.g20*map.g10 + map.det*map.g21*map.g11 + map.det*map.g22*map.g12;
241
 
    const real G0_2_2 = map.det*map.g20*map.g20 + map.det*map.g21*map.g21 + map.det*map.g22*map.g22;
242
 
 
243
 
    // Compute element tensor
244
 
    block[0] = 1.666666666666664e-01*G0_0_0 + 1.666666666666664e-01*G0_0_1 + 1.666666666666664e-01*G0_0_2 + 1.666666666666664e-01*G0_1_0 + 1.666666666666665e-01*G0_1_1 + 1.666666666666665e-01*G0_1_2 + 1.666666666666664e-01*G0_2_0 + 1.666666666666665e-01*G0_2_1 + 1.666666666666665e-01*G0_2_2;
245
 
    block[1] = -1.666666666666664e-01*G0_0_0 - 1.666666666666664e-01*G0_0_1 - 1.666666666666664e-01*G0_0_2;
246
 
    block[2] = -1.666666666666664e-01*G0_1_0 - 1.666666666666665e-01*G0_1_1 - 1.666666666666665e-01*G0_1_2;
247
 
    block[3] = -1.666666666666664e-01*G0_2_0 - 1.666666666666665e-01*G0_2_1 - 1.666666666666665e-01*G0_2_2;
248
 
    block[4] = -1.666666666666664e-01*G0_0_0 - 1.666666666666664e-01*G0_1_0 - 1.666666666666664e-01*G0_2_0;
249
 
    block[5] = 1.666666666666664e-01*G0_0_0;
250
 
    block[6] = 1.666666666666664e-01*G0_1_0;
251
 
    block[7] = 1.666666666666664e-01*G0_2_0;
252
 
    block[8] = -1.666666666666664e-01*G0_0_1 - 1.666666666666665e-01*G0_1_1 - 1.666666666666665e-01*G0_2_1;
253
 
    block[9] = 1.666666666666664e-01*G0_0_1;
254
 
    block[10] = 1.666666666666665e-01*G0_1_1;
255
 
    block[11] = 1.666666666666665e-01*G0_2_1;
256
 
    block[12] = -1.666666666666664e-01*G0_0_2 - 1.666666666666665e-01*G0_1_2 - 1.666666666666665e-01*G0_2_2;
257
 
    block[13] = 1.666666666666664e-01*G0_0_2;
258
 
    block[14] = 1.666666666666665e-01*G0_1_2;
259
 
    block[15] = 1.666666666666665e-01*G0_2_2;
260
 
  }
261
 
 
262
 
};
 
27
 
 
28
  class TestElement;
 
29
 
 
30
  class TrialElement;
 
31
 
 
32
  BilinearForm();
 
33
  
 
34
 
 
35
  bool interior_contribution() const;
 
36
 
 
37
  void eval(real block[], const AffineMap& map, real det) const;
 
38
 
 
39
  bool boundary_contribution() const;
 
40
 
 
41
  void eval(real block[], const AffineMap& map, real det, unsigned int facet) const;
 
42
 
 
43
  bool interior_boundary_contribution() const;
 
44
 
 
45
  void eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const;
 
46
 
 
47
};
 
48
 
 
49
class BilinearForm::TestElement : public dolfin::FiniteElement
 
50
{
 
51
public:
 
52
 
 
53
  TestElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)
 
54
  {
 
55
    // Element is scalar, don't need to initialize tensordims
 
56
 
 
57
    // Element is simple, don't need to initialize subelements
 
58
  }
 
59
 
 
60
  ~TestElement()
 
61
  {
 
62
    if ( tensordims ) delete [] tensordims;
 
63
    if ( subelements )
 
64
    {
 
65
      for (unsigned int i = 0; i < elementdim(); i++)
 
66
        delete subelements[i];
 
67
      delete [] subelements;
 
68
    }
 
69
  }
 
70
 
 
71
  inline unsigned int spacedim() const
 
72
  {
 
73
    return 4;
 
74
  }
 
75
 
 
76
  inline unsigned int shapedim() const
 
77
  {
 
78
    return 3;
 
79
  }
 
80
 
 
81
  inline unsigned int tensordim(unsigned int i) const
 
82
  {
 
83
    dolfin_error("Element is scalar.");
 
84
    return 0;
 
85
  }
 
86
 
 
87
  inline unsigned int elementdim() const
 
88
  {
 
89
    return 1;
 
90
  }
 
91
 
 
92
  inline unsigned int rank() const
 
93
  {
 
94
    return 0;
 
95
  }
 
96
 
 
97
  void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
 
98
  {
 
99
    nodes[0] = cell.entities(0)[0];
 
100
    nodes[1] = cell.entities(0)[1];
 
101
    nodes[2] = cell.entities(0)[2];
 
102
    nodes[3] = cell.entities(0)[3];
 
103
  }
 
104
 
 
105
  void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
 
106
  {
 
107
    points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
108
    points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
109
    points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
 
110
    points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
 
111
    components[0] = 0;
 
112
    components[1] = 0;
 
113
    components[2] = 0;
 
114
    components[3] = 0;
 
115
  }
 
116
 
 
117
  void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const
 
118
  {
 
119
    // FIXME: Temporary fix for Lagrange elements
 
120
    vertex_nodes[0] = vertex;
 
121
  }
 
122
 
 
123
  const FiniteElement& operator[] (unsigned int i) const
 
124
  {
 
125
    return *this;
 
126
  }
 
127
 
 
128
  FiniteElement& operator[] (unsigned int i)
 
129
  {
 
130
    return *this;
 
131
  }
 
132
 
 
133
  FiniteElementSpec spec() const
 
134
  {
 
135
    FiniteElementSpec s("Lagrange", "tetrahedron", 1);
 
136
    return s;
 
137
  }
 
138
  
 
139
private:
 
140
 
 
141
  unsigned int* tensordims;
 
142
  FiniteElement** subelements;
 
143
 
 
144
};
 
145
 
 
146
class BilinearForm::TrialElement : public dolfin::FiniteElement
 
147
{
 
148
public:
 
149
 
 
150
  TrialElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)
 
151
  {
 
152
    // Element is scalar, don't need to initialize tensordims
 
153
 
 
154
    // Element is simple, don't need to initialize subelements
 
155
  }
 
156
 
 
157
  ~TrialElement()
 
158
  {
 
159
    if ( tensordims ) delete [] tensordims;
 
160
    if ( subelements )
 
161
    {
 
162
      for (unsigned int i = 0; i < elementdim(); i++)
 
163
        delete subelements[i];
 
164
      delete [] subelements;
 
165
    }
 
166
  }
 
167
 
 
168
  inline unsigned int spacedim() const
 
169
  {
 
170
    return 4;
 
171
  }
 
172
 
 
173
  inline unsigned int shapedim() const
 
174
  {
 
175
    return 3;
 
176
  }
 
177
 
 
178
  inline unsigned int tensordim(unsigned int i) const
 
179
  {
 
180
    dolfin_error("Element is scalar.");
 
181
    return 0;
 
182
  }
 
183
 
 
184
  inline unsigned int elementdim() const
 
185
  {
 
186
    return 1;
 
187
  }
 
188
 
 
189
  inline unsigned int rank() const
 
190
  {
 
191
    return 0;
 
192
  }
 
193
 
 
194
  void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
 
195
  {
 
196
    nodes[0] = cell.entities(0)[0];
 
197
    nodes[1] = cell.entities(0)[1];
 
198
    nodes[2] = cell.entities(0)[2];
 
199
    nodes[3] = cell.entities(0)[3];
 
200
  }
 
201
 
 
202
  void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
 
203
  {
 
204
    points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
205
    points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
206
    points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
 
207
    points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
 
208
    components[0] = 0;
 
209
    components[1] = 0;
 
210
    components[2] = 0;
 
211
    components[3] = 0;
 
212
  }
 
213
 
 
214
  void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const
 
215
  {
 
216
    // FIXME: Temporary fix for Lagrange elements
 
217
    vertex_nodes[0] = vertex;
 
218
  }
 
219
 
 
220
  const FiniteElement& operator[] (unsigned int i) const
 
221
  {
 
222
    return *this;
 
223
  }
 
224
 
 
225
  FiniteElement& operator[] (unsigned int i)
 
226
  {
 
227
    return *this;
 
228
  }
 
229
 
 
230
  FiniteElementSpec spec() const
 
231
  {
 
232
    FiniteElementSpec s("Lagrange", "tetrahedron", 1);
 
233
    return s;
 
234
  }
 
235
  
 
236
private:
 
237
 
 
238
  unsigned int* tensordims;
 
239
  FiniteElement** subelements;
 
240
 
 
241
};
 
242
 
 
243
BilinearForm::BilinearForm() : dolfin::BilinearForm(0)
 
244
{
 
245
  // Create finite element for test space
 
246
  _test = new TestElement();
 
247
 
 
248
  // Create finite element for trial space
 
249
  _trial = new TrialElement();
 
250
}
 
251
 
 
252
// Contribution from the interior
 
253
bool BilinearForm::interior_contribution() const { return true; }
 
254
 
 
255
void BilinearForm::eval(real block[], const AffineMap& map, real det) const
 
256
{
 
257
  // Compute geometry tensors
 
258
  const real G0_0_0 = det*map.g00*map.g00 + det*map.g01*map.g01 + det*map.g02*map.g02;
 
259
  const real G0_0_1 = det*map.g00*map.g10 + det*map.g01*map.g11 + det*map.g02*map.g12;
 
260
  const real G0_0_2 = det*map.g00*map.g20 + det*map.g01*map.g21 + det*map.g02*map.g22;
 
261
  const real G0_1_0 = det*map.g10*map.g00 + det*map.g11*map.g01 + det*map.g12*map.g02;
 
262
  const real G0_1_1 = det*map.g10*map.g10 + det*map.g11*map.g11 + det*map.g12*map.g12;
 
263
  const real G0_1_2 = det*map.g10*map.g20 + det*map.g11*map.g21 + det*map.g12*map.g22;
 
264
  const real G0_2_0 = det*map.g20*map.g00 + det*map.g21*map.g01 + det*map.g22*map.g02;
 
265
  const real G0_2_1 = det*map.g20*map.g10 + det*map.g21*map.g11 + det*map.g22*map.g12;
 
266
  const real G0_2_2 = det*map.g20*map.g20 + det*map.g21*map.g21 + det*map.g22*map.g22;
 
267
 
 
268
  // Compute element tensor
 
269
  block[0] = 1.666666666666664e-01*G0_0_0 + 1.666666666666664e-01*G0_0_1 + 1.666666666666664e-01*G0_0_2 + 1.666666666666664e-01*G0_1_0 + 1.666666666666665e-01*G0_1_1 + 1.666666666666665e-01*G0_1_2 + 1.666666666666664e-01*G0_2_0 + 1.666666666666665e-01*G0_2_1 + 1.666666666666665e-01*G0_2_2;
 
270
  block[1] = -1.666666666666664e-01*G0_0_0 - 1.666666666666664e-01*G0_1_0 - 1.666666666666664e-01*G0_2_0;
 
271
  block[2] = -1.666666666666664e-01*G0_0_1 - 1.666666666666665e-01*G0_1_1 - 1.666666666666665e-01*G0_2_1;
 
272
  block[3] = -1.666666666666664e-01*G0_0_2 - 1.666666666666665e-01*G0_1_2 - 1.666666666666665e-01*G0_2_2;
 
273
  block[4] = -1.666666666666664e-01*G0_0_0 - 1.666666666666664e-01*G0_0_1 - 1.666666666666664e-01*G0_0_2;
 
274
  block[5] = 1.666666666666664e-01*G0_0_0;
 
275
  block[6] = 1.666666666666664e-01*G0_0_1;
 
276
  block[7] = 1.666666666666664e-01*G0_0_2;
 
277
  block[8] = -1.666666666666664e-01*G0_1_0 - 1.666666666666665e-01*G0_1_1 - 1.666666666666665e-01*G0_1_2;
 
278
  block[9] = 1.666666666666664e-01*G0_1_0;
 
279
  block[10] = 1.666666666666665e-01*G0_1_1;
 
280
  block[11] = 1.666666666666665e-01*G0_1_2;
 
281
  block[12] = -1.666666666666664e-01*G0_2_0 - 1.666666666666665e-01*G0_2_1 - 1.666666666666665e-01*G0_2_2;
 
282
  block[13] = 1.666666666666664e-01*G0_2_0;
 
283
  block[14] = 1.666666666666665e-01*G0_2_1;
 
284
  block[15] = 1.666666666666665e-01*G0_2_2;
 
285
}
 
286
 
 
287
// No contribution from the boundary
 
288
bool BilinearForm::boundary_contribution() const { return false; }
 
289
 
 
290
void BilinearForm::eval(real block[], const AffineMap& map, real det, unsigned int facet) const {}
 
291
 
 
292
// No contribution from interior boundaries
 
293
bool BilinearForm::interior_boundary_contribution() const { return false; }
 
294
 
 
295
void BilinearForm::eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const {}
263
296
 
264
297
/// This class contains the form to be evaluated, including
265
298
/// contributions from the interior and boundary of the domain.
267
300
class LinearForm : public dolfin::LinearForm
268
301
{
269
302
public:
270
 
  
271
 
  class TestElement : public dolfin::FiniteElement
272
 
  {
273
 
  public:
274
 
  
275
 
    TestElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)
276
 
    {
277
 
      // Element is scalar, don't need to initialize tensordims
278
 
  
279
 
      // Element is simple, don't need to initialize subelements
280
 
    }
281
 
  
282
 
    ~TestElement()
283
 
    {
284
 
      if ( tensordims ) delete [] tensordims;
285
 
      if ( subelements )
286
 
      {
287
 
        for (unsigned int i = 0; i < elementdim(); i++)
288
 
          delete subelements[i];
289
 
        delete [] subelements;
290
 
      }
291
 
    }
292
 
  
293
 
    inline unsigned int spacedim() const
294
 
    {
295
 
      return 4;
296
 
    }
297
 
  
298
 
    inline unsigned int shapedim() const
299
 
    {
300
 
      return 3;
301
 
    }
302
 
  
303
 
    inline unsigned int tensordim(unsigned int i) const
304
 
    {
305
 
      dolfin_error("Element is scalar.");
306
 
      return 0;
307
 
    }
308
 
  
309
 
    inline unsigned int elementdim() const
310
 
    {
311
 
      return 1;
312
 
    }
313
 
  
314
 
    inline unsigned int rank() const
315
 
    {
316
 
      return 0;
317
 
    }
318
 
  
319
 
    void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
320
 
    {
321
 
      nodes[0] = cell.vertexID(0);
322
 
      nodes[1] = cell.vertexID(1);
323
 
      nodes[2] = cell.vertexID(2);
324
 
      nodes[3] = cell.vertexID(3);
325
 
    }
326
 
  
327
 
    void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
328
 
    {
329
 
      points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
330
 
      points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
331
 
      points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
332
 
      points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
333
 
      components[0] = 0;
334
 
      components[1] = 0;
335
 
      components[2] = 0;
336
 
      components[3] = 0;
337
 
    }
338
 
  
339
 
    void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const
340
 
    {
341
 
      // FIXME: Temporary fix for Lagrange elements
342
 
      values[0] = x[vertex];
343
 
    }
344
 
  
345
 
    const FiniteElement& operator[] (unsigned int i) const
346
 
    {
347
 
      return *this;
348
 
    }
349
 
  
350
 
    FiniteElement& operator[] (unsigned int i)
351
 
    {
352
 
      return *this;
353
 
    }
354
 
  
355
 
    FiniteElementSpec spec() const
356
 
    {
357
 
      FiniteElementSpec s("Lagrange", "tetrahedron", 1);
358
 
      return s;
359
 
    }
360
 
    
361
 
  private:
362
 
  
363
 
    unsigned int* tensordims;
364
 
    FiniteElement** subelements;
365
 
  
366
 
  };
367
 
    
368
 
  class FunctionElement_0 : public dolfin::FiniteElement
369
 
  {
370
 
  public:
371
 
  
372
 
    FunctionElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)
373
 
    {
374
 
      // Element is scalar, don't need to initialize tensordims
375
 
  
376
 
      // Element is simple, don't need to initialize subelements
377
 
    }
378
 
  
379
 
    ~FunctionElement_0()
380
 
    {
381
 
      if ( tensordims ) delete [] tensordims;
382
 
      if ( subelements )
383
 
      {
384
 
        for (unsigned int i = 0; i < elementdim(); i++)
385
 
          delete subelements[i];
386
 
        delete [] subelements;
387
 
      }
388
 
    }
389
 
  
390
 
    inline unsigned int spacedim() const
391
 
    {
392
 
      return 4;
393
 
    }
394
 
  
395
 
    inline unsigned int shapedim() const
396
 
    {
397
 
      return 3;
398
 
    }
399
 
  
400
 
    inline unsigned int tensordim(unsigned int i) const
401
 
    {
402
 
      dolfin_error("Element is scalar.");
403
 
      return 0;
404
 
    }
405
 
  
406
 
    inline unsigned int elementdim() const
407
 
    {
408
 
      return 1;
409
 
    }
410
 
  
411
 
    inline unsigned int rank() const
412
 
    {
413
 
      return 0;
414
 
    }
415
 
  
416
 
    void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
417
 
    {
418
 
      nodes[0] = cell.vertexID(0);
419
 
      nodes[1] = cell.vertexID(1);
420
 
      nodes[2] = cell.vertexID(2);
421
 
      nodes[3] = cell.vertexID(3);
422
 
    }
423
 
  
424
 
    void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
425
 
    {
426
 
      points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
427
 
      points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
428
 
      points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
429
 
      points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
430
 
      components[0] = 0;
431
 
      components[1] = 0;
432
 
      components[2] = 0;
433
 
      components[3] = 0;
434
 
    }
435
 
  
436
 
    void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const
437
 
    {
438
 
      // FIXME: Temporary fix for Lagrange elements
439
 
      values[0] = x[vertex];
440
 
    }
441
 
  
442
 
    const FiniteElement& operator[] (unsigned int i) const
443
 
    {
444
 
      return *this;
445
 
    }
446
 
  
447
 
    FiniteElement& operator[] (unsigned int i)
448
 
    {
449
 
      return *this;
450
 
    }
451
 
  
452
 
    FiniteElementSpec spec() const
453
 
    {
454
 
      FiniteElementSpec s("Lagrange", "tetrahedron", 1);
455
 
      return s;
456
 
    }
457
 
    
458
 
  private:
459
 
  
460
 
    unsigned int* tensordims;
461
 
    FiniteElement** subelements;
462
 
  
463
 
  };
464
 
  
465
 
  LinearForm(Function& w0) : dolfin::LinearForm(1)
466
 
  {
467
 
    // Create finite element for test space
468
 
    _test = new TestElement();
469
 
 
470
 
    // Add functions
471
 
    add(w0, new FunctionElement_0());
472
 
  }
473
 
 
474
 
  void eval(real block[], const AffineMap& map) const
475
 
  {
476
 
    // Compute coefficients
477
 
    const real c0_0 = c[0][0];
478
 
    const real c0_1 = c[0][1];
479
 
    const real c0_2 = c[0][2];
480
 
    const real c0_3 = c[0][3];
481
 
 
482
 
    // Compute geometry tensors
483
 
    const real G0_0 = map.det*c0_0;
484
 
    const real G0_1 = map.det*c0_1;
485
 
    const real G0_2 = map.det*c0_2;
486
 
    const real G0_3 = map.det*c0_3;
487
 
 
488
 
    // Compute element tensor
489
 
    block[0] = 1.666666666666662e-02*G0_0 + 8.333333333333309e-03*G0_1 + 8.333333333333307e-03*G0_2 + 8.333333333333312e-03*G0_3;
490
 
    block[1] = 8.333333333333309e-03*G0_0 + 1.666666666666662e-02*G0_1 + 8.333333333333311e-03*G0_2 + 8.333333333333312e-03*G0_3;
491
 
    block[2] = 8.333333333333307e-03*G0_0 + 8.333333333333311e-03*G0_1 + 1.666666666666662e-02*G0_2 + 8.333333333333311e-03*G0_3;
492
 
    block[3] = 8.333333333333312e-03*G0_0 + 8.333333333333312e-03*G0_1 + 8.333333333333312e-03*G0_2 + 1.666666666666662e-02*G0_3;
493
 
  }
494
 
 
495
 
};
 
303
 
 
304
  class TestElement;
 
305
 
 
306
  class FunctionElement_0;
 
307
 
 
308
  LinearForm(Function& w0);
 
309
  
 
310
 
 
311
  bool interior_contribution() const;
 
312
 
 
313
  void eval(real block[], const AffineMap& map, real det) const;
 
314
 
 
315
  bool boundary_contribution() const;
 
316
 
 
317
  void eval(real block[], const AffineMap& map, real det, unsigned int facet) const;
 
318
 
 
319
  bool interior_boundary_contribution() const;
 
320
 
 
321
  void eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const;
 
322
 
 
323
};
 
324
 
 
325
class LinearForm::TestElement : public dolfin::FiniteElement
 
326
{
 
327
public:
 
328
 
 
329
  TestElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)
 
330
  {
 
331
    // Element is scalar, don't need to initialize tensordims
 
332
 
 
333
    // Element is simple, don't need to initialize subelements
 
334
  }
 
335
 
 
336
  ~TestElement()
 
337
  {
 
338
    if ( tensordims ) delete [] tensordims;
 
339
    if ( subelements )
 
340
    {
 
341
      for (unsigned int i = 0; i < elementdim(); i++)
 
342
        delete subelements[i];
 
343
      delete [] subelements;
 
344
    }
 
345
  }
 
346
 
 
347
  inline unsigned int spacedim() const
 
348
  {
 
349
    return 4;
 
350
  }
 
351
 
 
352
  inline unsigned int shapedim() const
 
353
  {
 
354
    return 3;
 
355
  }
 
356
 
 
357
  inline unsigned int tensordim(unsigned int i) const
 
358
  {
 
359
    dolfin_error("Element is scalar.");
 
360
    return 0;
 
361
  }
 
362
 
 
363
  inline unsigned int elementdim() const
 
364
  {
 
365
    return 1;
 
366
  }
 
367
 
 
368
  inline unsigned int rank() const
 
369
  {
 
370
    return 0;
 
371
  }
 
372
 
 
373
  void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
 
374
  {
 
375
    nodes[0] = cell.entities(0)[0];
 
376
    nodes[1] = cell.entities(0)[1];
 
377
    nodes[2] = cell.entities(0)[2];
 
378
    nodes[3] = cell.entities(0)[3];
 
379
  }
 
380
 
 
381
  void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
 
382
  {
 
383
    points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
384
    points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
385
    points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
 
386
    points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
 
387
    components[0] = 0;
 
388
    components[1] = 0;
 
389
    components[2] = 0;
 
390
    components[3] = 0;
 
391
  }
 
392
 
 
393
  void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const
 
394
  {
 
395
    // FIXME: Temporary fix for Lagrange elements
 
396
    vertex_nodes[0] = vertex;
 
397
  }
 
398
 
 
399
  const FiniteElement& operator[] (unsigned int i) const
 
400
  {
 
401
    return *this;
 
402
  }
 
403
 
 
404
  FiniteElement& operator[] (unsigned int i)
 
405
  {
 
406
    return *this;
 
407
  }
 
408
 
 
409
  FiniteElementSpec spec() const
 
410
  {
 
411
    FiniteElementSpec s("Lagrange", "tetrahedron", 1);
 
412
    return s;
 
413
  }
 
414
  
 
415
private:
 
416
 
 
417
  unsigned int* tensordims;
 
418
  FiniteElement** subelements;
 
419
 
 
420
};
 
421
 
 
422
class LinearForm::FunctionElement_0 : public dolfin::FiniteElement
 
423
{
 
424
public:
 
425
 
 
426
  FunctionElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)
 
427
  {
 
428
    // Element is scalar, don't need to initialize tensordims
 
429
 
 
430
    // Element is simple, don't need to initialize subelements
 
431
  }
 
432
 
 
433
  ~FunctionElement_0()
 
434
  {
 
435
    if ( tensordims ) delete [] tensordims;
 
436
    if ( subelements )
 
437
    {
 
438
      for (unsigned int i = 0; i < elementdim(); i++)
 
439
        delete subelements[i];
 
440
      delete [] subelements;
 
441
    }
 
442
  }
 
443
 
 
444
  inline unsigned int spacedim() const
 
445
  {
 
446
    return 4;
 
447
  }
 
448
 
 
449
  inline unsigned int shapedim() const
 
450
  {
 
451
    return 3;
 
452
  }
 
453
 
 
454
  inline unsigned int tensordim(unsigned int i) const
 
455
  {
 
456
    dolfin_error("Element is scalar.");
 
457
    return 0;
 
458
  }
 
459
 
 
460
  inline unsigned int elementdim() const
 
461
  {
 
462
    return 1;
 
463
  }
 
464
 
 
465
  inline unsigned int rank() const
 
466
  {
 
467
    return 0;
 
468
  }
 
469
 
 
470
  void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
 
471
  {
 
472
    nodes[0] = cell.entities(0)[0];
 
473
    nodes[1] = cell.entities(0)[1];
 
474
    nodes[2] = cell.entities(0)[2];
 
475
    nodes[3] = cell.entities(0)[3];
 
476
  }
 
477
 
 
478
  void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
 
479
  {
 
480
    points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
481
    points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
 
482
    points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
 
483
    points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
 
484
    components[0] = 0;
 
485
    components[1] = 0;
 
486
    components[2] = 0;
 
487
    components[3] = 0;
 
488
  }
 
489
 
 
490
  void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const
 
491
  {
 
492
    // FIXME: Temporary fix for Lagrange elements
 
493
    vertex_nodes[0] = vertex;
 
494
  }
 
495
 
 
496
  const FiniteElement& operator[] (unsigned int i) const
 
497
  {
 
498
    return *this;
 
499
  }
 
500
 
 
501
  FiniteElement& operator[] (unsigned int i)
 
502
  {
 
503
    return *this;
 
504
  }
 
505
 
 
506
  FiniteElementSpec spec() const
 
507
  {
 
508
    FiniteElementSpec s("Lagrange", "tetrahedron", 1);
 
509
    return s;
 
510
  }
 
511
  
 
512
private:
 
513
 
 
514
  unsigned int* tensordims;
 
515
  FiniteElement** subelements;
 
516
 
 
517
};
 
518
 
 
519
LinearForm::LinearForm(Function& w0) : dolfin::LinearForm(1)
 
520
{
 
521
  // Create finite element for test space
 
522
  _test = new TestElement();
 
523
 
 
524
  // Add functions
 
525
  initFunction(0, w0, new FunctionElement_0());
 
526
}
 
527
 
 
528
// Contribution from the interior
 
529
bool LinearForm::interior_contribution() const { return true; }
 
530
 
 
531
void LinearForm::eval(real block[], const AffineMap& map, real det) const
 
532
{
 
533
  // Compute coefficients
 
534
  const real c0_0 = c[0][0];
 
535
  const real c0_1 = c[0][1];
 
536
  const real c0_2 = c[0][2];
 
537
  const real c0_3 = c[0][3];
 
538
 
 
539
  // Compute geometry tensors
 
540
  const real G0_0 = det*c0_0;
 
541
  const real G0_1 = det*c0_1;
 
542
  const real G0_2 = det*c0_2;
 
543
  const real G0_3 = det*c0_3;
 
544
 
 
545
  // Compute element tensor
 
546
  block[0] = 1.666666666666662e-02*G0_0 + 8.333333333333309e-03*G0_1 + 8.333333333333307e-03*G0_2 + 8.333333333333312e-03*G0_3;
 
547
  block[1] = 8.333333333333309e-03*G0_0 + 1.666666666666662e-02*G0_1 + 8.333333333333311e-03*G0_2 + 8.333333333333312e-03*G0_3;
 
548
  block[2] = 8.333333333333307e-03*G0_0 + 8.333333333333311e-03*G0_1 + 1.666666666666662e-02*G0_2 + 8.333333333333312e-03*G0_3;
 
549
  block[3] = 8.333333333333312e-03*G0_0 + 8.333333333333312e-03*G0_1 + 8.333333333333311e-03*G0_2 + 1.666666666666662e-02*G0_3;
 
550
}
 
551
 
 
552
// No contribution from the boundary
 
553
bool LinearForm::boundary_contribution() const { return false; }
 
554
 
 
555
void LinearForm::eval(real block[], const AffineMap& map, real det, unsigned int facet) const {}
 
556
 
 
557
// No contribution from interior boundaries
 
558
bool LinearForm::interior_boundary_contribution() const { return false; }
 
559
 
 
560
void LinearForm::eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const {}
496
561
 
497
562
} }
498
563