← Back to branch summary

~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to bench/fem/assembly/cpp/forms/THStokes2D.h

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2014-09-22 14:35:34 UTC
  • mfrom: (1.1.17) (19.1.23 sid)
  • Revision ID: package-import@ubuntu.com-20140922143534-0yi89jyuqbgdxwm9
Tags: 1.4.0+dfsg-4
http://bugs.debian.org/760658
* debian/control: Disable libcgal-dev on i386, mipsel and sparc.
* debian/rules: Remove bad directives in pkg-config file dolfin.pc
  (closes: #760658).
* Remove debian/libdolfin-dev.lintian-overrides.

Show diffs side-by-side

added added

removed removed

Lines of Context:
bench/fem/assembly/cpp/forms/THStokes2D.h
1
 
// This code conforms with the UFC specification version 2.3.0
2
 
// and was automatically generated by FFC version 1.3.0.
 
1
// This code conforms with the UFC specification version 1.4.0
 
2
// and was automatically generated by FFC version 1.4.0.
3
3
//
4
4
// This code was generated with the option '-l dolfin' and
5
5
// contains DOLFIN-specific wrappers that depend on DOLFIN.
55
55
  /// Return a string identifying the finite element
56
56
  virtual const char* signature() const
57
57
  {
58
 
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
 
58
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 2, None)";
59
59
  }
60
60
 
61
61
  /// Return the cell shape
94
94
    return 1;
95
95
  }
96
96
 
97
 
  /// Evaluate basis function i at given point x in cell
 
97
  /// Evaluate basis function i at given point x in cell (actual implementation)
 
98
  static void _evaluate_basis(std::size_t i,
 
99
                              double* values,
 
100
                              const double* x,
 
101
                              const double* vertex_coordinates,
 
102
                              int cell_orientation)
 
103
  {
 
104
    // Compute Jacobian
 
105
    double J[4];
 
106
    compute_jacobian_triangle_2d(J, vertex_coordinates);
 
107
    
 
108
    // Compute Jacobian inverse and determinant
 
109
    double K[4];
 
110
    double detJ;
 
111
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
 
112
    
 
113
    
 
114
    // Compute constants
 
115
    const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
 
116
    const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
 
117
    
 
118
    // Get coordinates and map to the reference (FIAT) element
 
119
    double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
 
120
    double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
 
121
    
 
122
    // Reset values
 
123
    *values = 0.0;
 
124
    switch (i)
 
125
    {
 
126
    case 0:
 
127
      {
 
128
        
 
129
      // Array of basisvalues
 
130
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
131
      
 
132
      // Declare helper variables
 
133
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
134
      double tmp1 = (1.0 - Y)/2.0;
 
135
      double tmp2 = tmp1*tmp1;
 
136
      
 
137
      // Compute basisvalues
 
138
      basisvalues[0] = 1.0;
 
139
      basisvalues[1] = tmp0;
 
140
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
 
141
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
142
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
 
143
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
 
144
      basisvalues[0] *= std::sqrt(0.5);
 
145
      basisvalues[2] *= std::sqrt(1.0);
 
146
      basisvalues[5] *= std::sqrt(1.5);
 
147
      basisvalues[1] *= std::sqrt(3.0);
 
148
      basisvalues[4] *= std::sqrt(4.5);
 
149
      basisvalues[3] *= std::sqrt(7.5);
 
150
      
 
151
      // Table(s) of coefficients
 
152
      static const double coefficients0[6] = \
 
153
      {0.0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817};
 
154
      
 
155
      // Compute value(s)
 
156
      for (unsigned int r = 0; r < 6; r++)
 
157
      {
 
158
        *values += coefficients0[r]*basisvalues[r];
 
159
      } // end loop over 'r'
 
160
        break;
 
161
      }
 
162
    case 1:
 
163
      {
 
164
        
 
165
      // Array of basisvalues
 
166
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
167
      
 
168
      // Declare helper variables
 
169
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
170
      double tmp1 = (1.0 - Y)/2.0;
 
171
      double tmp2 = tmp1*tmp1;
 
172
      
 
173
      // Compute basisvalues
 
174
      basisvalues[0] = 1.0;
 
175
      basisvalues[1] = tmp0;
 
176
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
 
177
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
178
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
 
179
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
 
180
      basisvalues[0] *= std::sqrt(0.5);
 
181
      basisvalues[2] *= std::sqrt(1.0);
 
182
      basisvalues[5] *= std::sqrt(1.5);
 
183
      basisvalues[1] *= std::sqrt(3.0);
 
184
      basisvalues[4] *= std::sqrt(4.5);
 
185
      basisvalues[3] *= std::sqrt(7.5);
 
186
      
 
187
      // Table(s) of coefficients
 
188
      static const double coefficients0[6] = \
 
189
      {0.0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818};
 
190
      
 
191
      // Compute value(s)
 
192
      for (unsigned int r = 0; r < 6; r++)
 
193
      {
 
194
        *values += coefficients0[r]*basisvalues[r];
 
195
      } // end loop over 'r'
 
196
        break;
 
197
      }
 
198
    case 2:
 
199
      {
 
200
        
 
201
      // Array of basisvalues
 
202
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
203
      
 
204
      // Declare helper variables
 
205
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
206
      double tmp1 = (1.0 - Y)/2.0;
 
207
      double tmp2 = tmp1*tmp1;
 
208
      
 
209
      // Compute basisvalues
 
210
      basisvalues[0] = 1.0;
 
211
      basisvalues[1] = tmp0;
 
212
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
 
213
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
214
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
 
215
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
 
216
      basisvalues[0] *= std::sqrt(0.5);
 
217
      basisvalues[2] *= std::sqrt(1.0);
 
218
      basisvalues[5] *= std::sqrt(1.5);
 
219
      basisvalues[1] *= std::sqrt(3.0);
 
220
      basisvalues[4] *= std::sqrt(4.5);
 
221
      basisvalues[3] *= std::sqrt(7.5);
 
222
      
 
223
      // Table(s) of coefficients
 
224
      static const double coefficients0[6] = \
 
225
      {0.0, 0.0, 0.2, 0.0, 0.0, 0.163299316185545};
 
226
      
 
227
      // Compute value(s)
 
228
      for (unsigned int r = 0; r < 6; r++)
 
229
      {
 
230
        *values += coefficients0[r]*basisvalues[r];
 
231
      } // end loop over 'r'
 
232
        break;
 
233
      }
 
234
    case 3:
 
235
      {
 
236
        
 
237
      // Array of basisvalues
 
238
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
239
      
 
240
      // Declare helper variables
 
241
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
242
      double tmp1 = (1.0 - Y)/2.0;
 
243
      double tmp2 = tmp1*tmp1;
 
244
      
 
245
      // Compute basisvalues
 
246
      basisvalues[0] = 1.0;
 
247
      basisvalues[1] = tmp0;
 
248
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
 
249
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
250
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
 
251
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
 
252
      basisvalues[0] *= std::sqrt(0.5);
 
253
      basisvalues[2] *= std::sqrt(1.0);
 
254
      basisvalues[5] *= std::sqrt(1.5);
 
255
      basisvalues[1] *= std::sqrt(3.0);
 
256
      basisvalues[4] *= std::sqrt(4.5);
 
257
      basisvalues[3] *= std::sqrt(7.5);
 
258
      
 
259
      // Table(s) of coefficients
 
260
      static const double coefficients0[6] = \
 
261
      {0.471404520791032, 0.23094010767585, 0.133333333333333, 0.0, 0.188561808316413, -0.163299316185545};
 
262
      
 
263
      // Compute value(s)
 
264
      for (unsigned int r = 0; r < 6; r++)
 
265
      {
 
266
        *values += coefficients0[r]*basisvalues[r];
 
267
      } // end loop over 'r'
 
268
        break;
 
269
      }
 
270
    case 4:
 
271
      {
 
272
        
 
273
      // Array of basisvalues
 
274
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
275
      
 
276
      // Declare helper variables
 
277
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
278
      double tmp1 = (1.0 - Y)/2.0;
 
279
      double tmp2 = tmp1*tmp1;
 
280
      
 
281
      // Compute basisvalues
 
282
      basisvalues[0] = 1.0;
 
283
      basisvalues[1] = tmp0;
 
284
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
 
285
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
286
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
 
287
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
 
288
      basisvalues[0] *= std::sqrt(0.5);
 
289
      basisvalues[2] *= std::sqrt(1.0);
 
290
      basisvalues[5] *= std::sqrt(1.5);
 
291
      basisvalues[1] *= std::sqrt(3.0);
 
292
      basisvalues[4] *= std::sqrt(4.5);
 
293
      basisvalues[3] *= std::sqrt(7.5);
 
294
      
 
295
      // Table(s) of coefficients
 
296
      static const double coefficients0[6] = \
 
297
      {0.471404520791032, -0.23094010767585, 0.133333333333333, 0.0, -0.188561808316413, -0.163299316185545};
 
298
      
 
299
      // Compute value(s)
 
300
      for (unsigned int r = 0; r < 6; r++)
 
301
      {
 
302
        *values += coefficients0[r]*basisvalues[r];
 
303
      } // end loop over 'r'
 
304
        break;
 
305
      }
 
306
    case 5:
 
307
      {
 
308
        
 
309
      // Array of basisvalues
 
310
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
311
      
 
312
      // Declare helper variables
 
313
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
 
314
      double tmp1 = (1.0 - Y)/2.0;
 
315
      double tmp2 = tmp1*tmp1;
 
316
      
 
317
      // Compute basisvalues
 
318
      basisvalues[0] = 1.0;
 
319
      basisvalues[1] = tmp0;
 
320
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
 
321
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
 
322
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
 
323
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
 
324
      basisvalues[0] *= std::sqrt(0.5);
 
325
      basisvalues[2] *= std::sqrt(1.0);
 
326
      basisvalues[5] *= std::sqrt(1.5);
 
327
      basisvalues[1] *= std::sqrt(3.0);
 
328
      basisvalues[4] *= std::sqrt(4.5);
 
329
      basisvalues[3] *= std::sqrt(7.5);
 
330
      
 
331
      // Table(s) of coefficients
 
332
      static const double coefficients0[6] = \
 
333
      {0.471404520791032, 0.0, -0.266666666666667, -0.243432247780074, 0.0, 0.0544331053951817};
 
334
      
 
335
      // Compute value(s)
 
336
      for (unsigned int r = 0; r < 6; r++)
 
337
      {
 
338
        *values += coefficients0[r]*basisvalues[r];
 
339
      } // end loop over 'r'
 
340
        break;
 
341
      }
 
342
    }
 
343
    
 
344
  }
 
345
 
 
346
  /// Evaluate basis function i at given point x in cell (non-static member function)
98
347
  virtual void evaluate_basis(std::size_t i,
99
348
                              double* values,
100
349
                              const double* x,
101
350
                              const double* vertex_coordinates,
102
351
                              int cell_orientation) const
103
352
  {
104
 
    // Compute Jacobian
105
 
    double J[4];
106
 
    compute_jacobian_triangle_2d(J, vertex_coordinates);
107
 
    
108
 
    // Compute Jacobian inverse and determinant
109
 
    double K[4];
110
 
    double detJ;
111
 
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
112
 
    
113
 
    
114
 
    // Compute constants
115
 
    const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
116
 
    const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
117
 
    
118
 
    // Get coordinates and map to the reference (FIAT) element
119
 
    double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
120
 
    double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
121
 
    
122
 
    // Reset values
123
 
    *values = 0.0;
124
 
    switch (i)
125
 
    {
126
 
    case 0:
127
 
      {
128
 
        
129
 
      // Array of basisvalues
130
 
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
131
 
      
132
 
      // Declare helper variables
133
 
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
134
 
      double tmp1 = (1.0 - Y)/2.0;
135
 
      double tmp2 = tmp1*tmp1;
136
 
      
137
 
      // Compute basisvalues
138
 
      basisvalues[0] = 1.0;
139
 
      basisvalues[1] = tmp0;
140
 
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
141
 
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
142
 
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
143
 
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
144
 
      basisvalues[0] *= std::sqrt(0.5);
145
 
      basisvalues[2] *= std::sqrt(1.0);
146
 
      basisvalues[5] *= std::sqrt(1.5);
147
 
      basisvalues[1] *= std::sqrt(3.0);
148
 
      basisvalues[4] *= std::sqrt(4.5);
149
 
      basisvalues[3] *= std::sqrt(7.5);
150
 
      
151
 
      // Table(s) of coefficients
152
 
      static const double coefficients0[6] = \
153
 
      {0.0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817};
154
 
      
155
 
      // Compute value(s)
156
 
      for (unsigned int r = 0; r < 6; r++)
157
 
      {
158
 
        *values += coefficients0[r]*basisvalues[r];
159
 
      }// end loop over 'r'
160
 
        break;
161
 
      }
162
 
    case 1:
163
 
      {
164
 
        
165
 
      // Array of basisvalues
166
 
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
167
 
      
168
 
      // Declare helper variables
169
 
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
170
 
      double tmp1 = (1.0 - Y)/2.0;
171
 
      double tmp2 = tmp1*tmp1;
172
 
      
173
 
      // Compute basisvalues
174
 
      basisvalues[0] = 1.0;
175
 
      basisvalues[1] = tmp0;
176
 
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
177
 
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
178
 
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
179
 
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
180
 
      basisvalues[0] *= std::sqrt(0.5);
181
 
      basisvalues[2] *= std::sqrt(1.0);
182
 
      basisvalues[5] *= std::sqrt(1.5);
183
 
      basisvalues[1] *= std::sqrt(3.0);
184
 
      basisvalues[4] *= std::sqrt(4.5);
185
 
      basisvalues[3] *= std::sqrt(7.5);
186
 
      
187
 
      // Table(s) of coefficients
188
 
      static const double coefficients0[6] = \
189
 
      {0.0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818};
190
 
      
191
 
      // Compute value(s)
192
 
      for (unsigned int r = 0; r < 6; r++)
193
 
      {
194
 
        *values += coefficients0[r]*basisvalues[r];
195
 
      }// end loop over 'r'
196
 
        break;
197
 
      }
198
 
    case 2:
199
 
      {
200
 
        
201
 
      // Array of basisvalues
202
 
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
203
 
      
204
 
      // Declare helper variables
205
 
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
206
 
      double tmp1 = (1.0 - Y)/2.0;
207
 
      double tmp2 = tmp1*tmp1;
208
 
      
209
 
      // Compute basisvalues
210
 
      basisvalues[0] = 1.0;
211
 
      basisvalues[1] = tmp0;
212
 
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
213
 
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
214
 
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
215
 
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
216
 
      basisvalues[0] *= std::sqrt(0.5);
217
 
      basisvalues[2] *= std::sqrt(1.0);
218
 
      basisvalues[5] *= std::sqrt(1.5);
219
 
      basisvalues[1] *= std::sqrt(3.0);
220
 
      basisvalues[4] *= std::sqrt(4.5);
221
 
      basisvalues[3] *= std::sqrt(7.5);
222
 
      
223
 
      // Table(s) of coefficients
224
 
      static const double coefficients0[6] = \
225
 
      {0.0, 0.0, 0.2, 0.0, 0.0, 0.163299316185545};
226
 
      
227
 
      // Compute value(s)
228
 
      for (unsigned int r = 0; r < 6; r++)
229
 
      {
230
 
        *values += coefficients0[r]*basisvalues[r];
231
 
      }// end loop over 'r'
232
 
        break;
233
 
      }
234
 
    case 3:
235
 
      {
236
 
        
237
 
      // Array of basisvalues
238
 
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
239
 
      
240
 
      // Declare helper variables
241
 
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
242
 
      double tmp1 = (1.0 - Y)/2.0;
243
 
      double tmp2 = tmp1*tmp1;
244
 
      
245
 
      // Compute basisvalues
246
 
      basisvalues[0] = 1.0;
247
 
      basisvalues[1] = tmp0;
248
 
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
249
 
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
250
 
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
251
 
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
252
 
      basisvalues[0] *= std::sqrt(0.5);
253
 
      basisvalues[2] *= std::sqrt(1.0);
254
 
      basisvalues[5] *= std::sqrt(1.5);
255
 
      basisvalues[1] *= std::sqrt(3.0);
256
 
      basisvalues[4] *= std::sqrt(4.5);
257
 
      basisvalues[3] *= std::sqrt(7.5);
258
 
      
259
 
      // Table(s) of coefficients
260
 
      static const double coefficients0[6] = \
261
 
      {0.471404520791032, 0.23094010767585, 0.133333333333333, 0.0, 0.188561808316413, -0.163299316185545};
262
 
      
263
 
      // Compute value(s)
264
 
      for (unsigned int r = 0; r < 6; r++)
265
 
      {
266
 
        *values += coefficients0[r]*basisvalues[r];
267
 
      }// end loop over 'r'
268
 
        break;
269
 
      }
270
 
    case 4:
271
 
      {
272
 
        
273
 
      // Array of basisvalues
274
 
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
275
 
      
276
 
      // Declare helper variables
277
 
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
278
 
      double tmp1 = (1.0 - Y)/2.0;
279
 
      double tmp2 = tmp1*tmp1;
280
 
      
281
 
      // Compute basisvalues
282
 
      basisvalues[0] = 1.0;
283
 
      basisvalues[1] = tmp0;
284
 
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
285
 
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
286
 
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
287
 
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
288
 
      basisvalues[0] *= std::sqrt(0.5);
289
 
      basisvalues[2] *= std::sqrt(1.0);
290
 
      basisvalues[5] *= std::sqrt(1.5);
291
 
      basisvalues[1] *= std::sqrt(3.0);
292
 
      basisvalues[4] *= std::sqrt(4.5);
293
 
      basisvalues[3] *= std::sqrt(7.5);
294
 
      
295
 
      // Table(s) of coefficients
296
 
      static const double coefficients0[6] = \
297
 
      {0.471404520791032, -0.23094010767585, 0.133333333333333, 0.0, -0.188561808316413, -0.163299316185545};
298
 
      
299
 
      // Compute value(s)
300
 
      for (unsigned int r = 0; r < 6; r++)
301
 
      {
302
 
        *values += coefficients0[r]*basisvalues[r];
303
 
      }// end loop over 'r'
304
 
        break;
305
 
      }
306
 
    case 5:
307
 
      {
308
 
        
309
 
      // Array of basisvalues
310
 
      double basisvalues[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
311
 
      
312
 
      // Declare helper variables
313
 
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
314
 
      double tmp1 = (1.0 - Y)/2.0;
315
 
      double tmp2 = tmp1*tmp1;
316
 
      
317
 
      // Compute basisvalues
318
 
      basisvalues[0] = 1.0;
319
 
      basisvalues[1] = tmp0;
320
 
      basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
321
 
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
322
 
      basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
323
 
      basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
324
 
      basisvalues[0] *= std::sqrt(0.5);
325
 
      basisvalues[2] *= std::sqrt(1.0);
326
 
      basisvalues[5] *= std::sqrt(1.5);
327
 
      basisvalues[1] *= std::sqrt(3.0);
328
 
      basisvalues[4] *= std::sqrt(4.5);
329
 
      basisvalues[3] *= std::sqrt(7.5);
330
 
      
331
 
      // Table(s) of coefficients
332
 
      static const double coefficients0[6] = \
333
 
      {0.471404520791032, 0.0, -0.266666666666667, -0.243432247780074, 0.0, 0.0544331053951817};
334
 
      
335
 
      // Compute value(s)
336
 
      for (unsigned int r = 0; r < 6; r++)
337
 
      {
338
 
        *values += coefficients0[r]*basisvalues[r];
339
 
      }// end loop over 'r'
340
 
        break;
341
 
      }
342
 
    }
343
 
    
 
353
    _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
344
354
  }
345
355
 
346
 
  /// Evaluate all basis functions at given point x in cell
347
 
  virtual void evaluate_basis_all(double* values,
 
356
  /// Evaluate all basis functions at given point x in cell (actual implementation)
 
357
  static void _evaluate_basis_all(double* values,
348
358
                                  const double* x,
349
359
                                  const double* vertex_coordinates,
350
 
                                  int cell_orientation) const
 
360
                                  int cell_orientation)
351
361
  {
352
362
    // Helper variable to hold values of a single dof.
353
363
    double dof_values = 0.0;
355
365
    // Loop dofs and call evaluate_basis
356
366
    for (unsigned int r = 0; r < 6; r++)
357
367
    {
358
 
      evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
 
368
      _evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
359
369
      values[r] = dof_values;
360
 
    }// end loop over 'r'
361
 
  }
362
 
 
363
 
  /// Evaluate order n derivatives of basis function i at given point x in cell
364
 
  virtual void evaluate_basis_derivatives(std::size_t i,
 
370
    } // end loop over 'r'
 
371
  }
 
372
 
 
373
  /// Evaluate all basis functions at given point x in cell (non-static member function)
 
374
  virtual void evaluate_basis_all(double* values,
 
375
                                  const double* x,
 
376
                                  const double* vertex_coordinates,
 
377
                                  int cell_orientation) const
 
378
  {
 
379
    _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
380
  }
 
381
 
 
382
  /// Evaluate order n derivatives of basis function i at given point x in cell (actual implementation)
 
383
  static void _evaluate_basis_derivatives(std::size_t i,
365
384
                                          std::size_t n,
366
385
                                          double* values,
367
386
                                          const double* x,
368
387
                                          const double* vertex_coordinates,
369
 
                                          int cell_orientation) const
 
388
                                          int cell_orientation)
370
389
  {
371
390
    
372
391
    // Compute number of derivatives.
374
393
    for (unsigned int r = 0; r < n; r++)
375
394
    {
376
395
      num_derivatives *= 2;
377
 
    }// end loop over 'r'
 
396
    } // end loop over 'r'
378
397
    
379
398
    // Reset values. Assuming that values is always an array.
380
399
    for (unsigned int r = 0; r < num_derivatives; r++)
381
400
    {
382
401
      values[r] = 0.0;
383
 
    }// end loop over 'r'
 
402
    } // end loop over 'r'
384
403
    
385
404
    // Call evaluate_basis if order of derivatives is equal to zero.
386
405
    if (n == 0)
387
406
    {
388
 
      evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
 
407
      _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
389
408
      return ;
390
409
    }
391
410
    
514
533
      for (unsigned int r = 0; r < 4; r++)
515
534
      {
516
535
        derivatives[r] = 0.0;
517
 
      }// end loop over 'r'
 
536
      } // end loop over 'r'
518
537
      
519
538
      // Declare derivative matrix (of polynomial basis).
520
539
      double dmats[6][6] = \
548
567
            dmats[t][u] = 1.0;
549
568
            }
550
569
            
551
 
          }// end loop over 'u'
552
 
        }// end loop over 't'
 
570
          } // end loop over 'u'
 
571
        } // end loop over 't'
553
572
        
554
573
        // Looping derivative order to generate dmats.
555
574
        for (unsigned int s = 0; s < n; s++)
561
580
            {
562
581
              dmats_old[t][u] = dmats[t][u];
563
582
              dmats[t][u] = 0.0;
564
 
            }// end loop over 'u'
565
 
          }// end loop over 't'
 
583
            } // end loop over 'u'
 
584
          } // end loop over 't'
566
585
          
567
586
          // Update dmats using an inner product.
568
587
          if (combinations[r][s] == 0)
574
593
              for (unsigned int tu = 0; tu < 6; tu++)
575
594
              {
576
595
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
577
 
              }// end loop over 'tu'
578
 
            }// end loop over 'u'
579
 
          }// end loop over 't'
 
596
              } // end loop over 'tu'
 
597
            } // end loop over 'u'
 
598
          } // end loop over 't'
580
599
          }
581
600
          
582
601
          if (combinations[r][s] == 1)
588
607
              for (unsigned int tu = 0; tu < 6; tu++)
589
608
              {
590
609
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
591
 
              }// end loop over 'tu'
592
 
            }// end loop over 'u'
593
 
          }// end loop over 't'
 
610
              } // end loop over 'tu'
 
611
            } // end loop over 'u'
 
612
          } // end loop over 't'
594
613
          }
595
614
          
596
 
        }// end loop over 's'
 
615
        } // end loop over 's'
597
616
        for (unsigned int s = 0; s < 6; s++)
598
617
        {
599
618
          for (unsigned int t = 0; t < 6; t++)
600
619
          {
601
620
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
602
 
          }// end loop over 't'
603
 
        }// end loop over 's'
604
 
      }// end loop over 'r'
 
621
          } // end loop over 't'
 
622
        } // end loop over 's'
 
623
      } // end loop over 'r'
605
624
      
606
625
      // Transform derivatives back to physical element
607
626
      for (unsigned int r = 0; r < num_derivatives; r++)
609
628
        for (unsigned int s = 0; s < num_derivatives; s++)
610
629
        {
611
630
          values[r] += transform[r][s]*derivatives[s];
612
 
        }// end loop over 's'
613
 
      }// end loop over 'r'
 
631
        } // end loop over 's'
 
632
      } // end loop over 'r'
614
633
        break;
615
634
      }
616
635
    case 1:
665
684
      for (unsigned int r = 0; r < 4; r++)
666
685
      {
667
686
        derivatives[r] = 0.0;
668
 
      }// end loop over 'r'
 
687
      } // end loop over 'r'
669
688
      
670
689
      // Declare derivative matrix (of polynomial basis).
671
690
      double dmats[6][6] = \
699
718
            dmats[t][u] = 1.0;
700
719
            }
701
720
            
702
 
          }// end loop over 'u'
703
 
        }// end loop over 't'
 
721
          } // end loop over 'u'
 
722
        } // end loop over 't'
704
723
        
705
724
        // Looping derivative order to generate dmats.
706
725
        for (unsigned int s = 0; s < n; s++)
712
731
            {
713
732
              dmats_old[t][u] = dmats[t][u];
714
733
              dmats[t][u] = 0.0;
715
 
            }// end loop over 'u'
716
 
          }// end loop over 't'
 
734
            } // end loop over 'u'
 
735
          } // end loop over 't'
717
736
          
718
737
          // Update dmats using an inner product.
719
738
          if (combinations[r][s] == 0)
725
744
              for (unsigned int tu = 0; tu < 6; tu++)
726
745
              {
727
746
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
728
 
              }// end loop over 'tu'
729
 
            }// end loop over 'u'
730
 
          }// end loop over 't'
 
747
              } // end loop over 'tu'
 
748
            } // end loop over 'u'
 
749
          } // end loop over 't'
731
750
          }
732
751
          
733
752
          if (combinations[r][s] == 1)
739
758
              for (unsigned int tu = 0; tu < 6; tu++)
740
759
              {
741
760
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
742
 
              }// end loop over 'tu'
743
 
            }// end loop over 'u'
744
 
          }// end loop over 't'
 
761
              } // end loop over 'tu'
 
762
            } // end loop over 'u'
 
763
          } // end loop over 't'
745
764
          }
746
765
          
747
 
        }// end loop over 's'
 
766
        } // end loop over 's'
748
767
        for (unsigned int s = 0; s < 6; s++)
749
768
        {
750
769
          for (unsigned int t = 0; t < 6; t++)
751
770
          {
752
771
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
753
 
          }// end loop over 't'
754
 
        }// end loop over 's'
755
 
      }// end loop over 'r'
 
772
          } // end loop over 't'
 
773
        } // end loop over 's'
 
774
      } // end loop over 'r'
756
775
      
757
776
      // Transform derivatives back to physical element
758
777
      for (unsigned int r = 0; r < num_derivatives; r++)
760
779
        for (unsigned int s = 0; s < num_derivatives; s++)
761
780
        {
762
781
          values[r] += transform[r][s]*derivatives[s];
763
 
        }// end loop over 's'
764
 
      }// end loop over 'r'
 
782
        } // end loop over 's'
 
783
      } // end loop over 'r'
765
784
        break;
766
785
      }
767
786
    case 2:
816
835
      for (unsigned int r = 0; r < 4; r++)
817
836
      {
818
837
        derivatives[r] = 0.0;
819
 
      }// end loop over 'r'
 
838
      } // end loop over 'r'
820
839
      
821
840
      // Declare derivative matrix (of polynomial basis).
822
841
      double dmats[6][6] = \
850
869
            dmats[t][u] = 1.0;
851
870
            }
852
871
            
853
 
          }// end loop over 'u'
854
 
        }// end loop over 't'
 
872
          } // end loop over 'u'
 
873
        } // end loop over 't'
855
874
        
856
875
        // Looping derivative order to generate dmats.
857
876
        for (unsigned int s = 0; s < n; s++)
863
882
            {
864
883
              dmats_old[t][u] = dmats[t][u];
865
884
              dmats[t][u] = 0.0;
866
 
            }// end loop over 'u'
867
 
          }// end loop over 't'
 
885
            } // end loop over 'u'
 
886
          } // end loop over 't'
868
887
          
869
888
          // Update dmats using an inner product.
870
889
          if (combinations[r][s] == 0)
876
895
              for (unsigned int tu = 0; tu < 6; tu++)
877
896
              {
878
897
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
879
 
              }// end loop over 'tu'
880
 
            }// end loop over 'u'
881
 
          }// end loop over 't'
 
898
              } // end loop over 'tu'
 
899
            } // end loop over 'u'
 
900
          } // end loop over 't'
882
901
          }
883
902
          
884
903
          if (combinations[r][s] == 1)
890
909
              for (unsigned int tu = 0; tu < 6; tu++)
891
910
              {
892
911
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
893
 
              }// end loop over 'tu'
894
 
            }// end loop over 'u'
895
 
          }// end loop over 't'
 
912
              } // end loop over 'tu'
 
913
            } // end loop over 'u'
 
914
          } // end loop over 't'
896
915
          }
897
916
          
898
 
        }// end loop over 's'
 
917
        } // end loop over 's'
899
918
        for (unsigned int s = 0; s < 6; s++)
900
919
        {
901
920
          for (unsigned int t = 0; t < 6; t++)
902
921
          {
903
922
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
904
 
          }// end loop over 't'
905
 
        }// end loop over 's'
906
 
      }// end loop over 'r'
 
923
          } // end loop over 't'
 
924
        } // end loop over 's'
 
925
      } // end loop over 'r'
907
926
      
908
927
      // Transform derivatives back to physical element
909
928
      for (unsigned int r = 0; r < num_derivatives; r++)
911
930
        for (unsigned int s = 0; s < num_derivatives; s++)
912
931
        {
913
932
          values[r] += transform[r][s]*derivatives[s];
914
 
        }// end loop over 's'
915
 
      }// end loop over 'r'
 
933
        } // end loop over 's'
 
934
      } // end loop over 'r'
916
935
        break;
917
936
      }
918
937
    case 3:
967
986
      for (unsigned int r = 0; r < 4; r++)
968
987
      {
969
988
        derivatives[r] = 0.0;
970
 
      }// end loop over 'r'
 
989
      } // end loop over 'r'
971
990
      
972
991
      // Declare derivative matrix (of polynomial basis).
973
992
      double dmats[6][6] = \
1001
1020
            dmats[t][u] = 1.0;
1002
1021
            }
1003
1022
            
1004
 
          }// end loop over 'u'
1005
 
        }// end loop over 't'
 
1023
          } // end loop over 'u'
 
1024
        } // end loop over 't'
1006
1025
        
1007
1026
        // Looping derivative order to generate dmats.
1008
1027
        for (unsigned int s = 0; s < n; s++)
1014
1033
            {
1015
1034
              dmats_old[t][u] = dmats[t][u];
1016
1035
              dmats[t][u] = 0.0;
1017
 
            }// end loop over 'u'
1018
 
          }// end loop over 't'
 
1036
            } // end loop over 'u'
 
1037
          } // end loop over 't'
1019
1038
          
1020
1039
          // Update dmats using an inner product.
1021
1040
          if (combinations[r][s] == 0)
1027
1046
              for (unsigned int tu = 0; tu < 6; tu++)
1028
1047
              {
1029
1048
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1030
 
              }// end loop over 'tu'
1031
 
            }// end loop over 'u'
1032
 
          }// end loop over 't'
 
1049
              } // end loop over 'tu'
 
1050
            } // end loop over 'u'
 
1051
          } // end loop over 't'
1033
1052
          }
1034
1053
          
1035
1054
          if (combinations[r][s] == 1)
1041
1060
              for (unsigned int tu = 0; tu < 6; tu++)
1042
1061
              {
1043
1062
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1044
 
              }// end loop over 'tu'
1045
 
            }// end loop over 'u'
1046
 
          }// end loop over 't'
 
1063
              } // end loop over 'tu'
 
1064
            } // end loop over 'u'
 
1065
          } // end loop over 't'
1047
1066
          }
1048
1067
          
1049
 
        }// end loop over 's'
 
1068
        } // end loop over 's'
1050
1069
        for (unsigned int s = 0; s < 6; s++)
1051
1070
        {
1052
1071
          for (unsigned int t = 0; t < 6; t++)
1053
1072
          {
1054
1073
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1055
 
          }// end loop over 't'
1056
 
        }// end loop over 's'
1057
 
      }// end loop over 'r'
 
1074
          } // end loop over 't'
 
1075
        } // end loop over 's'
 
1076
      } // end loop over 'r'
1058
1077
      
1059
1078
      // Transform derivatives back to physical element
1060
1079
      for (unsigned int r = 0; r < num_derivatives; r++)
1062
1081
        for (unsigned int s = 0; s < num_derivatives; s++)
1063
1082
        {
1064
1083
          values[r] += transform[r][s]*derivatives[s];
1065
 
        }// end loop over 's'
1066
 
      }// end loop over 'r'
 
1084
        } // end loop over 's'
 
1085
      } // end loop over 'r'
1067
1086
        break;
1068
1087
      }
1069
1088
    case 4:
1118
1137
      for (unsigned int r = 0; r < 4; r++)
1119
1138
      {
1120
1139
        derivatives[r] = 0.0;
1121
 
      }// end loop over 'r'
 
1140
      } // end loop over 'r'
1122
1141
      
1123
1142
      // Declare derivative matrix (of polynomial basis).
1124
1143
      double dmats[6][6] = \
1152
1171
            dmats[t][u] = 1.0;
1153
1172
            }
1154
1173
            
1155
 
          }// end loop over 'u'
1156
 
        }// end loop over 't'
 
1174
          } // end loop over 'u'
 
1175
        } // end loop over 't'
1157
1176
        
1158
1177
        // Looping derivative order to generate dmats.
1159
1178
        for (unsigned int s = 0; s < n; s++)
1165
1184
            {
1166
1185
              dmats_old[t][u] = dmats[t][u];
1167
1186
              dmats[t][u] = 0.0;
1168
 
            }// end loop over 'u'
1169
 
          }// end loop over 't'
 
1187
            } // end loop over 'u'
 
1188
          } // end loop over 't'
1170
1189
          
1171
1190
          // Update dmats using an inner product.
1172
1191
          if (combinations[r][s] == 0)
1178
1197
              for (unsigned int tu = 0; tu < 6; tu++)
1179
1198
              {
1180
1199
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1181
 
              }// end loop over 'tu'
1182
 
            }// end loop over 'u'
1183
 
          }// end loop over 't'
 
1200
              } // end loop over 'tu'
 
1201
            } // end loop over 'u'
 
1202
          } // end loop over 't'
1184
1203
          }
1185
1204
          
1186
1205
          if (combinations[r][s] == 1)
1192
1211
              for (unsigned int tu = 0; tu < 6; tu++)
1193
1212
              {
1194
1213
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1195
 
              }// end loop over 'tu'
1196
 
            }// end loop over 'u'
1197
 
          }// end loop over 't'
 
1214
              } // end loop over 'tu'
 
1215
            } // end loop over 'u'
 
1216
          } // end loop over 't'
1198
1217
          }
1199
1218
          
1200
 
        }// end loop over 's'
 
1219
        } // end loop over 's'
1201
1220
        for (unsigned int s = 0; s < 6; s++)
1202
1221
        {
1203
1222
          for (unsigned int t = 0; t < 6; t++)
1204
1223
          {
1205
1224
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1206
 
          }// end loop over 't'
1207
 
        }// end loop over 's'
1208
 
      }// end loop over 'r'
 
1225
          } // end loop over 't'
 
1226
        } // end loop over 's'
 
1227
      } // end loop over 'r'
1209
1228
      
1210
1229
      // Transform derivatives back to physical element
1211
1230
      for (unsigned int r = 0; r < num_derivatives; r++)
1213
1232
        for (unsigned int s = 0; s < num_derivatives; s++)
1214
1233
        {
1215
1234
          values[r] += transform[r][s]*derivatives[s];
1216
 
        }// end loop over 's'
1217
 
      }// end loop over 'r'
 
1235
        } // end loop over 's'
 
1236
      } // end loop over 'r'
1218
1237
        break;
1219
1238
      }
1220
1239
    case 5:
1269
1288
      for (unsigned int r = 0; r < 4; r++)
1270
1289
      {
1271
1290
        derivatives[r] = 0.0;
1272
 
      }// end loop over 'r'
 
1291
      } // end loop over 'r'
1273
1292
      
1274
1293
      // Declare derivative matrix (of polynomial basis).
1275
1294
      double dmats[6][6] = \
1303
1322
            dmats[t][u] = 1.0;
1304
1323
            }
1305
1324
            
1306
 
          }// end loop over 'u'
1307
 
        }// end loop over 't'
 
1325
          } // end loop over 'u'
 
1326
        } // end loop over 't'
1308
1327
        
1309
1328
        // Looping derivative order to generate dmats.
1310
1329
        for (unsigned int s = 0; s < n; s++)
1316
1335
            {
1317
1336
              dmats_old[t][u] = dmats[t][u];
1318
1337
              dmats[t][u] = 0.0;
1319
 
            }// end loop over 'u'
1320
 
          }// end loop over 't'
 
1338
            } // end loop over 'u'
 
1339
          } // end loop over 't'
1321
1340
          
1322
1341
          // Update dmats using an inner product.
1323
1342
          if (combinations[r][s] == 0)
1329
1348
              for (unsigned int tu = 0; tu < 6; tu++)
1330
1349
              {
1331
1350
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1332
 
              }// end loop over 'tu'
1333
 
            }// end loop over 'u'
1334
 
          }// end loop over 't'
 
1351
              } // end loop over 'tu'
 
1352
            } // end loop over 'u'
 
1353
          } // end loop over 't'
1335
1354
          }
1336
1355
          
1337
1356
          if (combinations[r][s] == 1)
1343
1362
              for (unsigned int tu = 0; tu < 6; tu++)
1344
1363
              {
1345
1364
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1346
 
              }// end loop over 'tu'
1347
 
            }// end loop over 'u'
1348
 
          }// end loop over 't'
 
1365
              } // end loop over 'tu'
 
1366
            } // end loop over 'u'
 
1367
          } // end loop over 't'
1349
1368
          }
1350
1369
          
1351
 
        }// end loop over 's'
 
1370
        } // end loop over 's'
1352
1371
        for (unsigned int s = 0; s < 6; s++)
1353
1372
        {
1354
1373
          for (unsigned int t = 0; t < 6; t++)
1355
1374
          {
1356
1375
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1357
 
          }// end loop over 't'
1358
 
        }// end loop over 's'
1359
 
      }// end loop over 'r'
 
1376
          } // end loop over 't'
 
1377
        } // end loop over 's'
 
1378
      } // end loop over 'r'
1360
1379
      
1361
1380
      // Transform derivatives back to physical element
1362
1381
      for (unsigned int r = 0; r < num_derivatives; r++)
1364
1383
        for (unsigned int s = 0; s < num_derivatives; s++)
1365
1384
        {
1366
1385
          values[r] += transform[r][s]*derivatives[s];
1367
 
        }// end loop over 's'
1368
 
      }// end loop over 'r'
 
1386
        } // end loop over 's'
 
1387
      } // end loop over 'r'
1369
1388
        break;
1370
1389
      }
1371
1390
    }
1372
1391
    
1373
1392
  }
1374
1393
 
1375
 
  /// Evaluate order n derivatives of all basis functions at given point x in cell
1376
 
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
1394
  /// Evaluate order n derivatives of basis function i at given point x in cell (non-static member function)
 
1395
  virtual void evaluate_basis_derivatives(std::size_t i,
 
1396
                                          std::size_t n,
 
1397
                                          double* values,
 
1398
                                          const double* x,
 
1399
                                          const double* vertex_coordinates,
 
1400
                                          int cell_orientation) const
 
1401
  {
 
1402
    _evaluate_basis_derivatives(i, n, values, x, vertex_coordinates, cell_orientation);
 
1403
  }
 
1404
 
 
1405
  /// Evaluate order n derivatives of all basis functions at given point x in cell (actual implementation)
 
1406
  static void _evaluate_basis_derivatives_all(std::size_t n,
1377
1407
                                              double* values,
1378
1408
                                              const double* x,
1379
1409
                                              const double* vertex_coordinates,
1380
 
                                              int cell_orientation) const
 
1410
                                              int cell_orientation)
1381
1411
  {
1382
1412
    // Call evaluate_basis_all if order of derivatives is equal to zero.
1383
1413
    if (n == 0)
1384
1414
    {
1385
 
      evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
1415
      _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
1386
1416
      return ;
1387
1417
    }
1388
1418
    
1391
1421
    for (unsigned int r = 0; r < n; r++)
1392
1422
    {
1393
1423
      num_derivatives *= 2;
1394
 
    }// end loop over 'r'
 
1424
    } // end loop over 'r'
1395
1425
    
1396
1426
    // Set values equal to zero.
1397
1427
    for (unsigned int r = 0; r < 6; r++)
1399
1429
      for (unsigned int s = 0; s < num_derivatives; s++)
1400
1430
      {
1401
1431
        values[r*num_derivatives + s] = 0.0;
1402
 
      }// end loop over 's'
1403
 
    }// end loop over 'r'
 
1432
      } // end loop over 's'
 
1433
    } // end loop over 'r'
1404
1434
    
1405
1435
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
1406
1436
    if (n > 2)
1413
1443
    for (unsigned int r = 0; r < 4; r++)
1414
1444
    {
1415
1445
      dof_values[r] = 0.0;
1416
 
    }// end loop over 'r'
 
1446
    } // end loop over 'r'
1417
1447
    
1418
1448
    // Loop dofs and call evaluate_basis_derivatives.
1419
1449
    for (unsigned int r = 0; r < 6; r++)
1420
1450
    {
1421
 
      evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
 
1451
      _evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1422
1452
      for (unsigned int s = 0; s < num_derivatives; s++)
1423
1453
      {
1424
1454
        values[r*num_derivatives + s] = dof_values[s];
1425
 
      }// end loop over 's'
1426
 
    }// end loop over 'r'
 
1455
      } // end loop over 's'
 
1456
    } // end loop over 'r'
 
1457
  }
 
1458
 
 
1459
  /// Evaluate order n derivatives of all basis functions at given point x in cell (non-static member function)
 
1460
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
1461
                                              double* values,
 
1462
                                              const double* x,
 
1463
                                              const double* vertex_coordinates,
 
1464
                                              int cell_orientation) const
 
1465
  {
 
1466
    _evaluate_basis_derivatives_all(n, values, x, vertex_coordinates, cell_orientation);
1427
1467
  }
1428
1468
 
1429
1469
  /// Evaluate linear functional for dof i on the function f
1601
1641
  /// Return a string identifying the finite element
1602
1642
  virtual const char* signature() const
1603
1643
  {
1604
 
    return "VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None)";
 
1644
    return "VectorElement('Lagrange', Domain(Cell('triangle', 2)), 2, 2, None)";
1605
1645
  }
1606
1646
 
1607
1647
  /// Return the cell shape
1649
1689
    return 0;
1650
1690
  }
1651
1691
 
1652
 
  /// Evaluate basis function i at given point x in cell
1653
 
  virtual void evaluate_basis(std::size_t i,
 
1692
  /// Evaluate basis function i at given point x in cell (actual implementation)
 
1693
  static void _evaluate_basis(std::size_t i,
1654
1694
                              double* values,
1655
1695
                              const double* x,
1656
1696
                              const double* vertex_coordinates,
1657
 
                              int cell_orientation) const
 
1697
                              int cell_orientation)
1658
1698
  {
1659
1699
    // Compute Jacobian
1660
1700
    double J[4];
1712
1752
      for (unsigned int r = 0; r < 6; r++)
1713
1753
      {
1714
1754
        values[0] += coefficients0[r]*basisvalues[r];
1715
 
      }// end loop over 'r'
 
1755
      } // end loop over 'r'
1716
1756
        break;
1717
1757
      }
1718
1758
    case 1:
1748
1788
      for (unsigned int r = 0; r < 6; r++)
1749
1789
      {
1750
1790
        values[0] += coefficients0[r]*basisvalues[r];
1751
 
      }// end loop over 'r'
 
1791
      } // end loop over 'r'
1752
1792
        break;
1753
1793
      }
1754
1794
    case 2:
1784
1824
      for (unsigned int r = 0; r < 6; r++)
1785
1825
      {
1786
1826
        values[0] += coefficients0[r]*basisvalues[r];
1787
 
      }// end loop over 'r'
 
1827
      } // end loop over 'r'
1788
1828
        break;
1789
1829
      }
1790
1830
    case 3:
1820
1860
      for (unsigned int r = 0; r < 6; r++)
1821
1861
      {
1822
1862
        values[0] += coefficients0[r]*basisvalues[r];
1823
 
      }// end loop over 'r'
 
1863
      } // end loop over 'r'
1824
1864
        break;
1825
1865
      }
1826
1866
    case 4:
1856
1896
      for (unsigned int r = 0; r < 6; r++)
1857
1897
      {
1858
1898
        values[0] += coefficients0[r]*basisvalues[r];
1859
 
      }// end loop over 'r'
 
1899
      } // end loop over 'r'
1860
1900
        break;
1861
1901
      }
1862
1902
    case 5:
1892
1932
      for (unsigned int r = 0; r < 6; r++)
1893
1933
      {
1894
1934
        values[0] += coefficients0[r]*basisvalues[r];
1895
 
      }// end loop over 'r'
 
1935
      } // end loop over 'r'
1896
1936
        break;
1897
1937
      }
1898
1938
    case 6:
1928
1968
      for (unsigned int r = 0; r < 6; r++)
1929
1969
      {
1930
1970
        values[1] += coefficients0[r]*basisvalues[r];
1931
 
      }// end loop over 'r'
 
1971
      } // end loop over 'r'
1932
1972
        break;
1933
1973
      }
1934
1974
    case 7:
1964
2004
      for (unsigned int r = 0; r < 6; r++)
1965
2005
      {
1966
2006
        values[1] += coefficients0[r]*basisvalues[r];
1967
 
      }// end loop over 'r'
 
2007
      } // end loop over 'r'
1968
2008
        break;
1969
2009
      }
1970
2010
    case 8:
2000
2040
      for (unsigned int r = 0; r < 6; r++)
2001
2041
      {
2002
2042
        values[1] += coefficients0[r]*basisvalues[r];
2003
 
      }// end loop over 'r'
 
2043
      } // end loop over 'r'
2004
2044
        break;
2005
2045
      }
2006
2046
    case 9:
2036
2076
      for (unsigned int r = 0; r < 6; r++)
2037
2077
      {
2038
2078
        values[1] += coefficients0[r]*basisvalues[r];
2039
 
      }// end loop over 'r'
 
2079
      } // end loop over 'r'
2040
2080
        break;
2041
2081
      }
2042
2082
    case 10:
2072
2112
      for (unsigned int r = 0; r < 6; r++)
2073
2113
      {
2074
2114
        values[1] += coefficients0[r]*basisvalues[r];
2075
 
      }// end loop over 'r'
 
2115
      } // end loop over 'r'
2076
2116
        break;
2077
2117
      }
2078
2118
    case 11:
2108
2148
      for (unsigned int r = 0; r < 6; r++)
2109
2149
      {
2110
2150
        values[1] += coefficients0[r]*basisvalues[r];
2111
 
      }// end loop over 'r'
 
2151
      } // end loop over 'r'
2112
2152
        break;
2113
2153
      }
2114
2154
    }
2115
2155
    
2116
2156
  }
2117
2157
 
2118
 
  /// Evaluate all basis functions at given point x in cell
2119
 
  virtual void evaluate_basis_all(double* values,
 
2158
  /// Evaluate basis function i at given point x in cell (non-static member function)
 
2159
  virtual void evaluate_basis(std::size_t i,
 
2160
                              double* values,
 
2161
                              const double* x,
 
2162
                              const double* vertex_coordinates,
 
2163
                              int cell_orientation) const
 
2164
  {
 
2165
    _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
 
2166
  }
 
2167
 
 
2168
  /// Evaluate all basis functions at given point x in cell (actual implementation)
 
2169
  static void _evaluate_basis_all(double* values,
2120
2170
                                  const double* x,
2121
2171
                                  const double* vertex_coordinates,
2122
 
                                  int cell_orientation) const
 
2172
                                  int cell_orientation)
2123
2173
  {
2124
2174
    // Helper variable to hold values of a single dof.
2125
2175
    double dof_values[2] = {0.0, 0.0};
2127
2177
    // Loop dofs and call evaluate_basis
2128
2178
    for (unsigned int r = 0; r < 12; r++)
2129
2179
    {
2130
 
      evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
 
2180
      _evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
2131
2181
      for (unsigned int s = 0; s < 2; s++)
2132
2182
      {
2133
2183
        values[r*2 + s] = dof_values[s];
2134
 
      }// end loop over 's'
2135
 
    }// end loop over 'r'
2136
 
  }
2137
 
 
2138
 
  /// Evaluate order n derivatives of basis function i at given point x in cell
2139
 
  virtual void evaluate_basis_derivatives(std::size_t i,
 
2184
      } // end loop over 's'
 
2185
    } // end loop over 'r'
 
2186
  }
 
2187
 
 
2188
  /// Evaluate all basis functions at given point x in cell (non-static member function)
 
2189
  virtual void evaluate_basis_all(double* values,
 
2190
                                  const double* x,
 
2191
                                  const double* vertex_coordinates,
 
2192
                                  int cell_orientation) const
 
2193
  {
 
2194
    _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
2195
  }
 
2196
 
 
2197
  /// Evaluate order n derivatives of basis function i at given point x in cell (actual implementation)
 
2198
  static void _evaluate_basis_derivatives(std::size_t i,
2140
2199
                                          std::size_t n,
2141
2200
                                          double* values,
2142
2201
                                          const double* x,
2143
2202
                                          const double* vertex_coordinates,
2144
 
                                          int cell_orientation) const
 
2203
                                          int cell_orientation)
2145
2204
  {
2146
2205
    
2147
2206
    // Compute number of derivatives.
2149
2208
    for (unsigned int r = 0; r < n; r++)
2150
2209
    {
2151
2210
      num_derivatives *= 2;
2152
 
    }// end loop over 'r'
 
2211
    } // end loop over 'r'
2153
2212
    
2154
2213
    // Reset values. Assuming that values is always an array.
2155
2214
    for (unsigned int r = 0; r < 2*num_derivatives; r++)
2156
2215
    {
2157
2216
      values[r] = 0.0;
2158
 
    }// end loop over 'r'
 
2217
    } // end loop over 'r'
2159
2218
    
2160
2219
    // Call evaluate_basis if order of derivatives is equal to zero.
2161
2220
    if (n == 0)
2162
2221
    {
2163
 
      evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
 
2222
      _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
2164
2223
      return ;
2165
2224
    }
2166
2225
    
2289
2348
      for (unsigned int r = 0; r < 4; r++)
2290
2349
      {
2291
2350
        derivatives[r] = 0.0;
2292
 
      }// end loop over 'r'
 
2351
      } // end loop over 'r'
2293
2352
      
2294
2353
      // Declare derivative matrix (of polynomial basis).
2295
2354
      double dmats[6][6] = \
2323
2382
            dmats[t][u] = 1.0;
2324
2383
            }
2325
2384
            
2326
 
          }// end loop over 'u'
2327
 
        }// end loop over 't'
 
2385
          } // end loop over 'u'
 
2386
        } // end loop over 't'
2328
2387
        
2329
2388
        // Looping derivative order to generate dmats.
2330
2389
        for (unsigned int s = 0; s < n; s++)
2336
2395
            {
2337
2396
              dmats_old[t][u] = dmats[t][u];
2338
2397
              dmats[t][u] = 0.0;
2339
 
            }// end loop over 'u'
2340
 
          }// end loop over 't'
 
2398
            } // end loop over 'u'
 
2399
          } // end loop over 't'
2341
2400
          
2342
2401
          // Update dmats using an inner product.
2343
2402
          if (combinations[r][s] == 0)
2349
2408
              for (unsigned int tu = 0; tu < 6; tu++)
2350
2409
              {
2351
2410
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2352
 
              }// end loop over 'tu'
2353
 
            }// end loop over 'u'
2354
 
          }// end loop over 't'
 
2411
              } // end loop over 'tu'
 
2412
            } // end loop over 'u'
 
2413
          } // end loop over 't'
2355
2414
          }
2356
2415
          
2357
2416
          if (combinations[r][s] == 1)
2363
2422
              for (unsigned int tu = 0; tu < 6; tu++)
2364
2423
              {
2365
2424
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2366
 
              }// end loop over 'tu'
2367
 
            }// end loop over 'u'
2368
 
          }// end loop over 't'
 
2425
              } // end loop over 'tu'
 
2426
            } // end loop over 'u'
 
2427
          } // end loop over 't'
2369
2428
          }
2370
2429
          
2371
 
        }// end loop over 's'
 
2430
        } // end loop over 's'
2372
2431
        for (unsigned int s = 0; s < 6; s++)
2373
2432
        {
2374
2433
          for (unsigned int t = 0; t < 6; t++)
2375
2434
          {
2376
2435
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2377
 
          }// end loop over 't'
2378
 
        }// end loop over 's'
2379
 
      }// end loop over 'r'
 
2436
          } // end loop over 't'
 
2437
        } // end loop over 's'
 
2438
      } // end loop over 'r'
2380
2439
      
2381
2440
      // Transform derivatives back to physical element
2382
2441
      for (unsigned int r = 0; r < num_derivatives; r++)
2384
2443
        for (unsigned int s = 0; s < num_derivatives; s++)
2385
2444
        {
2386
2445
          values[r] += transform[r][s]*derivatives[s];
2387
 
        }// end loop over 's'
2388
 
      }// end loop over 'r'
 
2446
        } // end loop over 's'
 
2447
      } // end loop over 'r'
2389
2448
        break;
2390
2449
      }
2391
2450
    case 1:
2440
2499
      for (unsigned int r = 0; r < 4; r++)
2441
2500
      {
2442
2501
        derivatives[r] = 0.0;
2443
 
      }// end loop over 'r'
 
2502
      } // end loop over 'r'
2444
2503
      
2445
2504
      // Declare derivative matrix (of polynomial basis).
2446
2505
      double dmats[6][6] = \
2474
2533
            dmats[t][u] = 1.0;
2475
2534
            }
2476
2535
            
2477
 
          }// end loop over 'u'
2478
 
        }// end loop over 't'
 
2536
          } // end loop over 'u'
 
2537
        } // end loop over 't'
2479
2538
        
2480
2539
        // Looping derivative order to generate dmats.
2481
2540
        for (unsigned int s = 0; s < n; s++)
2487
2546
            {
2488
2547
              dmats_old[t][u] = dmats[t][u];
2489
2548
              dmats[t][u] = 0.0;
2490
 
            }// end loop over 'u'
2491
 
          }// end loop over 't'
 
2549
            } // end loop over 'u'
 
2550
          } // end loop over 't'
2492
2551
          
2493
2552
          // Update dmats using an inner product.
2494
2553
          if (combinations[r][s] == 0)
2500
2559
              for (unsigned int tu = 0; tu < 6; tu++)
2501
2560
              {
2502
2561
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2503
 
              }// end loop over 'tu'
2504
 
            }// end loop over 'u'
2505
 
          }// end loop over 't'
 
2562
              } // end loop over 'tu'
 
2563
            } // end loop over 'u'
 
2564
          } // end loop over 't'
2506
2565
          }
2507
2566
          
2508
2567
          if (combinations[r][s] == 1)
2514
2573
              for (unsigned int tu = 0; tu < 6; tu++)
2515
2574
              {
2516
2575
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2517
 
              }// end loop over 'tu'
2518
 
            }// end loop over 'u'
2519
 
          }// end loop over 't'
 
2576
              } // end loop over 'tu'
 
2577
            } // end loop over 'u'
 
2578
          } // end loop over 't'
2520
2579
          }
2521
2580
          
2522
 
        }// end loop over 's'
 
2581
        } // end loop over 's'
2523
2582
        for (unsigned int s = 0; s < 6; s++)
2524
2583
        {
2525
2584
          for (unsigned int t = 0; t < 6; t++)
2526
2585
          {
2527
2586
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2528
 
          }// end loop over 't'
2529
 
        }// end loop over 's'
2530
 
      }// end loop over 'r'
 
2587
          } // end loop over 't'
 
2588
        } // end loop over 's'
 
2589
      } // end loop over 'r'
2531
2590
      
2532
2591
      // Transform derivatives back to physical element
2533
2592
      for (unsigned int r = 0; r < num_derivatives; r++)
2535
2594
        for (unsigned int s = 0; s < num_derivatives; s++)
2536
2595
        {
2537
2596
          values[r] += transform[r][s]*derivatives[s];
2538
 
        }// end loop over 's'
2539
 
      }// end loop over 'r'
 
2597
        } // end loop over 's'
 
2598
      } // end loop over 'r'
2540
2599
        break;
2541
2600
      }
2542
2601
    case 2:
2591
2650
      for (unsigned int r = 0; r < 4; r++)
2592
2651
      {
2593
2652
        derivatives[r] = 0.0;
2594
 
      }// end loop over 'r'
 
2653
      } // end loop over 'r'
2595
2654
      
2596
2655
      // Declare derivative matrix (of polynomial basis).
2597
2656
      double dmats[6][6] = \
2625
2684
            dmats[t][u] = 1.0;
2626
2685
            }
2627
2686
            
2628
 
          }// end loop over 'u'
2629
 
        }// end loop over 't'
 
2687
          } // end loop over 'u'
 
2688
        } // end loop over 't'
2630
2689
        
2631
2690
        // Looping derivative order to generate dmats.
2632
2691
        for (unsigned int s = 0; s < n; s++)
2638
2697
            {
2639
2698
              dmats_old[t][u] = dmats[t][u];
2640
2699
              dmats[t][u] = 0.0;
2641
 
            }// end loop over 'u'
2642
 
          }// end loop over 't'
 
2700
            } // end loop over 'u'
 
2701
          } // end loop over 't'
2643
2702
          
2644
2703
          // Update dmats using an inner product.
2645
2704
          if (combinations[r][s] == 0)
2651
2710
              for (unsigned int tu = 0; tu < 6; tu++)
2652
2711
              {
2653
2712
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2654
 
              }// end loop over 'tu'
2655
 
            }// end loop over 'u'
2656
 
          }// end loop over 't'
 
2713
              } // end loop over 'tu'
 
2714
            } // end loop over 'u'
 
2715
          } // end loop over 't'
2657
2716
          }
2658
2717
          
2659
2718
          if (combinations[r][s] == 1)
2665
2724
              for (unsigned int tu = 0; tu < 6; tu++)
2666
2725
              {
2667
2726
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2668
 
              }// end loop over 'tu'
2669
 
            }// end loop over 'u'
2670
 
          }// end loop over 't'
 
2727
              } // end loop over 'tu'
 
2728
            } // end loop over 'u'
 
2729
          } // end loop over 't'
2671
2730
          }
2672
2731
          
2673
 
        }// end loop over 's'
 
2732
        } // end loop over 's'
2674
2733
        for (unsigned int s = 0; s < 6; s++)
2675
2734
        {
2676
2735
          for (unsigned int t = 0; t < 6; t++)
2677
2736
          {
2678
2737
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2679
 
          }// end loop over 't'
2680
 
        }// end loop over 's'
2681
 
      }// end loop over 'r'
 
2738
          } // end loop over 't'
 
2739
        } // end loop over 's'
 
2740
      } // end loop over 'r'
2682
2741
      
2683
2742
      // Transform derivatives back to physical element
2684
2743
      for (unsigned int r = 0; r < num_derivatives; r++)
2686
2745
        for (unsigned int s = 0; s < num_derivatives; s++)
2687
2746
        {
2688
2747
          values[r] += transform[r][s]*derivatives[s];
2689
 
        }// end loop over 's'
2690
 
      }// end loop over 'r'
 
2748
        } // end loop over 's'
 
2749
      } // end loop over 'r'
2691
2750
        break;
2692
2751
      }
2693
2752
    case 3:
2742
2801
      for (unsigned int r = 0; r < 4; r++)
2743
2802
      {
2744
2803
        derivatives[r] = 0.0;
2745
 
      }// end loop over 'r'
 
2804
      } // end loop over 'r'
2746
2805
      
2747
2806
      // Declare derivative matrix (of polynomial basis).
2748
2807
      double dmats[6][6] = \
2776
2835
            dmats[t][u] = 1.0;
2777
2836
            }
2778
2837
            
2779
 
          }// end loop over 'u'
2780
 
        }// end loop over 't'
 
2838
          } // end loop over 'u'
 
2839
        } // end loop over 't'
2781
2840
        
2782
2841
        // Looping derivative order to generate dmats.
2783
2842
        for (unsigned int s = 0; s < n; s++)
2789
2848
            {
2790
2849
              dmats_old[t][u] = dmats[t][u];
2791
2850
              dmats[t][u] = 0.0;
2792
 
            }// end loop over 'u'
2793
 
          }// end loop over 't'
 
2851
            } // end loop over 'u'
 
2852
          } // end loop over 't'
2794
2853
          
2795
2854
          // Update dmats using an inner product.
2796
2855
          if (combinations[r][s] == 0)
2802
2861
              for (unsigned int tu = 0; tu < 6; tu++)
2803
2862
              {
2804
2863
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2805
 
              }// end loop over 'tu'
2806
 
            }// end loop over 'u'
2807
 
          }// end loop over 't'
 
2864
              } // end loop over 'tu'
 
2865
            } // end loop over 'u'
 
2866
          } // end loop over 't'
2808
2867
          }
2809
2868
          
2810
2869
          if (combinations[r][s] == 1)
2816
2875
              for (unsigned int tu = 0; tu < 6; tu++)
2817
2876
              {
2818
2877
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2819
 
              }// end loop over 'tu'
2820
 
            }// end loop over 'u'
2821
 
          }// end loop over 't'
 
2878
              } // end loop over 'tu'
 
2879
            } // end loop over 'u'
 
2880
          } // end loop over 't'
2822
2881
          }
2823
2882
          
2824
 
        }// end loop over 's'
 
2883
        } // end loop over 's'
2825
2884
        for (unsigned int s = 0; s < 6; s++)
2826
2885
        {
2827
2886
          for (unsigned int t = 0; t < 6; t++)
2828
2887
          {
2829
2888
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2830
 
          }// end loop over 't'
2831
 
        }// end loop over 's'
2832
 
      }// end loop over 'r'
 
2889
          } // end loop over 't'
 
2890
        } // end loop over 's'
 
2891
      } // end loop over 'r'
2833
2892
      
2834
2893
      // Transform derivatives back to physical element
2835
2894
      for (unsigned int r = 0; r < num_derivatives; r++)
2837
2896
        for (unsigned int s = 0; s < num_derivatives; s++)
2838
2897
        {
2839
2898
          values[r] += transform[r][s]*derivatives[s];
2840
 
        }// end loop over 's'
2841
 
      }// end loop over 'r'
 
2899
        } // end loop over 's'
 
2900
      } // end loop over 'r'
2842
2901
        break;
2843
2902
      }
2844
2903
    case 4:
2893
2952
      for (unsigned int r = 0; r < 4; r++)
2894
2953
      {
2895
2954
        derivatives[r] = 0.0;
2896
 
      }// end loop over 'r'
 
2955
      } // end loop over 'r'
2897
2956
      
2898
2957
      // Declare derivative matrix (of polynomial basis).
2899
2958
      double dmats[6][6] = \
2927
2986
            dmats[t][u] = 1.0;
2928
2987
            }
2929
2988
            
2930
 
          }// end loop over 'u'
2931
 
        }// end loop over 't'
 
2989
          } // end loop over 'u'
 
2990
        } // end loop over 't'
2932
2991
        
2933
2992
        // Looping derivative order to generate dmats.
2934
2993
        for (unsigned int s = 0; s < n; s++)
2940
2999
            {
2941
3000
              dmats_old[t][u] = dmats[t][u];
2942
3001
              dmats[t][u] = 0.0;
2943
 
            }// end loop over 'u'
2944
 
          }// end loop over 't'
 
3002
            } // end loop over 'u'
 
3003
          } // end loop over 't'
2945
3004
          
2946
3005
          // Update dmats using an inner product.
2947
3006
          if (combinations[r][s] == 0)
2953
3012
              for (unsigned int tu = 0; tu < 6; tu++)
2954
3013
              {
2955
3014
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2956
 
              }// end loop over 'tu'
2957
 
            }// end loop over 'u'
2958
 
          }// end loop over 't'
 
3015
              } // end loop over 'tu'
 
3016
            } // end loop over 'u'
 
3017
          } // end loop over 't'
2959
3018
          }
2960
3019
          
2961
3020
          if (combinations[r][s] == 1)
2967
3026
              for (unsigned int tu = 0; tu < 6; tu++)
2968
3027
              {
2969
3028
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2970
 
              }// end loop over 'tu'
2971
 
            }// end loop over 'u'
2972
 
          }// end loop over 't'
 
3029
              } // end loop over 'tu'
 
3030
            } // end loop over 'u'
 
3031
          } // end loop over 't'
2973
3032
          }
2974
3033
          
2975
 
        }// end loop over 's'
 
3034
        } // end loop over 's'
2976
3035
        for (unsigned int s = 0; s < 6; s++)
2977
3036
        {
2978
3037
          for (unsigned int t = 0; t < 6; t++)
2979
3038
          {
2980
3039
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2981
 
          }// end loop over 't'
2982
 
        }// end loop over 's'
2983
 
      }// end loop over 'r'
 
3040
          } // end loop over 't'
 
3041
        } // end loop over 's'
 
3042
      } // end loop over 'r'
2984
3043
      
2985
3044
      // Transform derivatives back to physical element
2986
3045
      for (unsigned int r = 0; r < num_derivatives; r++)
2988
3047
        for (unsigned int s = 0; s < num_derivatives; s++)
2989
3048
        {
2990
3049
          values[r] += transform[r][s]*derivatives[s];
2991
 
        }// end loop over 's'
2992
 
      }// end loop over 'r'
 
3050
        } // end loop over 's'
 
3051
      } // end loop over 'r'
2993
3052
        break;
2994
3053
      }
2995
3054
    case 5:
3044
3103
      for (unsigned int r = 0; r < 4; r++)
3045
3104
      {
3046
3105
        derivatives[r] = 0.0;
3047
 
      }// end loop over 'r'
 
3106
      } // end loop over 'r'
3048
3107
      
3049
3108
      // Declare derivative matrix (of polynomial basis).
3050
3109
      double dmats[6][6] = \
3078
3137
            dmats[t][u] = 1.0;
3079
3138
            }
3080
3139
            
3081
 
          }// end loop over 'u'
3082
 
        }// end loop over 't'
 
3140
          } // end loop over 'u'
 
3141
        } // end loop over 't'
3083
3142
        
3084
3143
        // Looping derivative order to generate dmats.
3085
3144
        for (unsigned int s = 0; s < n; s++)
3091
3150
            {
3092
3151
              dmats_old[t][u] = dmats[t][u];
3093
3152
              dmats[t][u] = 0.0;
3094
 
            }// end loop over 'u'
3095
 
          }// end loop over 't'
 
3153
            } // end loop over 'u'
 
3154
          } // end loop over 't'
3096
3155
          
3097
3156
          // Update dmats using an inner product.
3098
3157
          if (combinations[r][s] == 0)
3104
3163
              for (unsigned int tu = 0; tu < 6; tu++)
3105
3164
              {
3106
3165
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3107
 
              }// end loop over 'tu'
3108
 
            }// end loop over 'u'
3109
 
          }// end loop over 't'
 
3166
              } // end loop over 'tu'
 
3167
            } // end loop over 'u'
 
3168
          } // end loop over 't'
3110
3169
          }
3111
3170
          
3112
3171
          if (combinations[r][s] == 1)
3118
3177
              for (unsigned int tu = 0; tu < 6; tu++)
3119
3178
              {
3120
3179
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3121
 
              }// end loop over 'tu'
3122
 
            }// end loop over 'u'
3123
 
          }// end loop over 't'
 
3180
              } // end loop over 'tu'
 
3181
            } // end loop over 'u'
 
3182
          } // end loop over 't'
3124
3183
          }
3125
3184
          
3126
 
        }// end loop over 's'
 
3185
        } // end loop over 's'
3127
3186
        for (unsigned int s = 0; s < 6; s++)
3128
3187
        {
3129
3188
          for (unsigned int t = 0; t < 6; t++)
3130
3189
          {
3131
3190
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3132
 
          }// end loop over 't'
3133
 
        }// end loop over 's'
3134
 
      }// end loop over 'r'
 
3191
          } // end loop over 't'
 
3192
        } // end loop over 's'
 
3193
      } // end loop over 'r'
3135
3194
      
3136
3195
      // Transform derivatives back to physical element
3137
3196
      for (unsigned int r = 0; r < num_derivatives; r++)
3139
3198
        for (unsigned int s = 0; s < num_derivatives; s++)
3140
3199
        {
3141
3200
          values[r] += transform[r][s]*derivatives[s];
3142
 
        }// end loop over 's'
3143
 
      }// end loop over 'r'
 
3201
        } // end loop over 's'
 
3202
      } // end loop over 'r'
3144
3203
        break;
3145
3204
      }
3146
3205
    case 6:
3195
3254
      for (unsigned int r = 0; r < 4; r++)
3196
3255
      {
3197
3256
        derivatives[r] = 0.0;
3198
 
      }// end loop over 'r'
 
3257
      } // end loop over 'r'
3199
3258
      
3200
3259
      // Declare derivative matrix (of polynomial basis).
3201
3260
      double dmats[6][6] = \
3229
3288
            dmats[t][u] = 1.0;
3230
3289
            }
3231
3290
            
3232
 
          }// end loop over 'u'
3233
 
        }// end loop over 't'
 
3291
          } // end loop over 'u'
 
3292
        } // end loop over 't'
3234
3293
        
3235
3294
        // Looping derivative order to generate dmats.
3236
3295
        for (unsigned int s = 0; s < n; s++)
3242
3301
            {
3243
3302
              dmats_old[t][u] = dmats[t][u];
3244
3303
              dmats[t][u] = 0.0;
3245
 
            }// end loop over 'u'
3246
 
          }// end loop over 't'
 
3304
            } // end loop over 'u'
 
3305
          } // end loop over 't'
3247
3306
          
3248
3307
          // Update dmats using an inner product.
3249
3308
          if (combinations[r][s] == 0)
3255
3314
              for (unsigned int tu = 0; tu < 6; tu++)
3256
3315
              {
3257
3316
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3258
 
              }// end loop over 'tu'
3259
 
            }// end loop over 'u'
3260
 
          }// end loop over 't'
 
3317
              } // end loop over 'tu'
 
3318
            } // end loop over 'u'
 
3319
          } // end loop over 't'
3261
3320
          }
3262
3321
          
3263
3322
          if (combinations[r][s] == 1)
3269
3328
              for (unsigned int tu = 0; tu < 6; tu++)
3270
3329
              {
3271
3330
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3272
 
              }// end loop over 'tu'
3273
 
            }// end loop over 'u'
3274
 
          }// end loop over 't'
 
3331
              } // end loop over 'tu'
 
3332
            } // end loop over 'u'
 
3333
          } // end loop over 't'
3275
3334
          }
3276
3335
          
3277
 
        }// end loop over 's'
 
3336
        } // end loop over 's'
3278
3337
        for (unsigned int s = 0; s < 6; s++)
3279
3338
        {
3280
3339
          for (unsigned int t = 0; t < 6; t++)
3281
3340
          {
3282
3341
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3283
 
          }// end loop over 't'
3284
 
        }// end loop over 's'
3285
 
      }// end loop over 'r'
 
3342
          } // end loop over 't'
 
3343
        } // end loop over 's'
 
3344
      } // end loop over 'r'
3286
3345
      
3287
3346
      // Transform derivatives back to physical element
3288
3347
      for (unsigned int r = 0; r < num_derivatives; r++)
3290
3349
        for (unsigned int s = 0; s < num_derivatives; s++)
3291
3350
        {
3292
3351
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
3293
 
        }// end loop over 's'
3294
 
      }// end loop over 'r'
 
3352
        } // end loop over 's'
 
3353
      } // end loop over 'r'
3295
3354
        break;
3296
3355
      }
3297
3356
    case 7:
3346
3405
      for (unsigned int r = 0; r < 4; r++)
3347
3406
      {
3348
3407
        derivatives[r] = 0.0;
3349
 
      }// end loop over 'r'
 
3408
      } // end loop over 'r'
3350
3409
      
3351
3410
      // Declare derivative matrix (of polynomial basis).
3352
3411
      double dmats[6][6] = \
3380
3439
            dmats[t][u] = 1.0;
3381
3440
            }
3382
3441
            
3383
 
          }// end loop over 'u'
3384
 
        }// end loop over 't'
 
3442
          } // end loop over 'u'
 
3443
        } // end loop over 't'
3385
3444
        
3386
3445
        // Looping derivative order to generate dmats.
3387
3446
        for (unsigned int s = 0; s < n; s++)
3393
3452
            {
3394
3453
              dmats_old[t][u] = dmats[t][u];
3395
3454
              dmats[t][u] = 0.0;
3396
 
            }// end loop over 'u'
3397
 
          }// end loop over 't'
 
3455
            } // end loop over 'u'
 
3456
          } // end loop over 't'
3398
3457
          
3399
3458
          // Update dmats using an inner product.
3400
3459
          if (combinations[r][s] == 0)
3406
3465
              for (unsigned int tu = 0; tu < 6; tu++)
3407
3466
              {
3408
3467
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3409
 
              }// end loop over 'tu'
3410
 
            }// end loop over 'u'
3411
 
          }// end loop over 't'
 
3468
              } // end loop over 'tu'
 
3469
            } // end loop over 'u'
 
3470
          } // end loop over 't'
3412
3471
          }
3413
3472
          
3414
3473
          if (combinations[r][s] == 1)
3420
3479
              for (unsigned int tu = 0; tu < 6; tu++)
3421
3480
              {
3422
3481
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3423
 
              }// end loop over 'tu'
3424
 
            }// end loop over 'u'
3425
 
          }// end loop over 't'
 
3482
              } // end loop over 'tu'
 
3483
            } // end loop over 'u'
 
3484
          } // end loop over 't'
3426
3485
          }
3427
3486
          
3428
 
        }// end loop over 's'
 
3487
        } // end loop over 's'
3429
3488
        for (unsigned int s = 0; s < 6; s++)
3430
3489
        {
3431
3490
          for (unsigned int t = 0; t < 6; t++)
3432
3491
          {
3433
3492
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3434
 
          }// end loop over 't'
3435
 
        }// end loop over 's'
3436
 
      }// end loop over 'r'
 
3493
          } // end loop over 't'
 
3494
        } // end loop over 's'
 
3495
      } // end loop over 'r'
3437
3496
      
3438
3497
      // Transform derivatives back to physical element
3439
3498
      for (unsigned int r = 0; r < num_derivatives; r++)
3441
3500
        for (unsigned int s = 0; s < num_derivatives; s++)
3442
3501
        {
3443
3502
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
3444
 
        }// end loop over 's'
3445
 
      }// end loop over 'r'
 
3503
        } // end loop over 's'
 
3504
      } // end loop over 'r'
3446
3505
        break;
3447
3506
      }
3448
3507
    case 8:
3497
3556
      for (unsigned int r = 0; r < 4; r++)
3498
3557
      {
3499
3558
        derivatives[r] = 0.0;
3500
 
      }// end loop over 'r'
 
3559
      } // end loop over 'r'
3501
3560
      
3502
3561
      // Declare derivative matrix (of polynomial basis).
3503
3562
      double dmats[6][6] = \
3531
3590
            dmats[t][u] = 1.0;
3532
3591
            }
3533
3592
            
3534
 
          }// end loop over 'u'
3535
 
        }// end loop over 't'
 
3593
          } // end loop over 'u'
 
3594
        } // end loop over 't'
3536
3595
        
3537
3596
        // Looping derivative order to generate dmats.
3538
3597
        for (unsigned int s = 0; s < n; s++)
3544
3603
            {
3545
3604
              dmats_old[t][u] = dmats[t][u];
3546
3605
              dmats[t][u] = 0.0;
3547
 
            }// end loop over 'u'
3548
 
          }// end loop over 't'
 
3606
            } // end loop over 'u'
 
3607
          } // end loop over 't'
3549
3608
          
3550
3609
          // Update dmats using an inner product.
3551
3610
          if (combinations[r][s] == 0)
3557
3616
              for (unsigned int tu = 0; tu < 6; tu++)
3558
3617
              {
3559
3618
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3560
 
              }// end loop over 'tu'
3561
 
            }// end loop over 'u'
3562
 
          }// end loop over 't'
 
3619
              } // end loop over 'tu'
 
3620
            } // end loop over 'u'
 
3621
          } // end loop over 't'
3563
3622
          }
3564
3623
          
3565
3624
          if (combinations[r][s] == 1)
3571
3630
              for (unsigned int tu = 0; tu < 6; tu++)
3572
3631
              {
3573
3632
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3574
 
              }// end loop over 'tu'
3575
 
            }// end loop over 'u'
3576
 
          }// end loop over 't'
 
3633
              } // end loop over 'tu'
 
3634
            } // end loop over 'u'
 
3635
          } // end loop over 't'
3577
3636
          }
3578
3637
          
3579
 
        }// end loop over 's'
 
3638
        } // end loop over 's'
3580
3639
        for (unsigned int s = 0; s < 6; s++)
3581
3640
        {
3582
3641
          for (unsigned int t = 0; t < 6; t++)
3583
3642
          {
3584
3643
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3585
 
          }// end loop over 't'
3586
 
        }// end loop over 's'
3587
 
      }// end loop over 'r'
 
3644
          } // end loop over 't'
 
3645
        } // end loop over 's'
 
3646
      } // end loop over 'r'
3588
3647
      
3589
3648
      // Transform derivatives back to physical element
3590
3649
      for (unsigned int r = 0; r < num_derivatives; r++)
3592
3651
        for (unsigned int s = 0; s < num_derivatives; s++)
3593
3652
        {
3594
3653
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
3595
 
        }// end loop over 's'
3596
 
      }// end loop over 'r'
 
3654
        } // end loop over 's'
 
3655
      } // end loop over 'r'
3597
3656
        break;
3598
3657
      }
3599
3658
    case 9:
3648
3707
      for (unsigned int r = 0; r < 4; r++)
3649
3708
      {
3650
3709
        derivatives[r] = 0.0;
3651
 
      }// end loop over 'r'
 
3710
      } // end loop over 'r'
3652
3711
      
3653
3712
      // Declare derivative matrix (of polynomial basis).
3654
3713
      double dmats[6][6] = \
3682
3741
            dmats[t][u] = 1.0;
3683
3742
            }
3684
3743
            
3685
 
          }// end loop over 'u'
3686
 
        }// end loop over 't'
 
3744
          } // end loop over 'u'
 
3745
        } // end loop over 't'
3687
3746
        
3688
3747
        // Looping derivative order to generate dmats.
3689
3748
        for (unsigned int s = 0; s < n; s++)
3695
3754
            {
3696
3755
              dmats_old[t][u] = dmats[t][u];
3697
3756
              dmats[t][u] = 0.0;
3698
 
            }// end loop over 'u'
3699
 
          }// end loop over 't'
 
3757
            } // end loop over 'u'
 
3758
          } // end loop over 't'
3700
3759
          
3701
3760
          // Update dmats using an inner product.
3702
3761
          if (combinations[r][s] == 0)
3708
3767
              for (unsigned int tu = 0; tu < 6; tu++)
3709
3768
              {
3710
3769
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3711
 
              }// end loop over 'tu'
3712
 
            }// end loop over 'u'
3713
 
          }// end loop over 't'
 
3770
              } // end loop over 'tu'
 
3771
            } // end loop over 'u'
 
3772
          } // end loop over 't'
3714
3773
          }
3715
3774
          
3716
3775
          if (combinations[r][s] == 1)
3722
3781
              for (unsigned int tu = 0; tu < 6; tu++)
3723
3782
              {
3724
3783
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3725
 
              }// end loop over 'tu'
3726
 
            }// end loop over 'u'
3727
 
          }// end loop over 't'
 
3784
              } // end loop over 'tu'
 
3785
            } // end loop over 'u'
 
3786
          } // end loop over 't'
3728
3787
          }
3729
3788
          
3730
 
        }// end loop over 's'
 
3789
        } // end loop over 's'
3731
3790
        for (unsigned int s = 0; s < 6; s++)
3732
3791
        {
3733
3792
          for (unsigned int t = 0; t < 6; t++)
3734
3793
          {
3735
3794
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3736
 
          }// end loop over 't'
3737
 
        }// end loop over 's'
3738
 
      }// end loop over 'r'
 
3795
          } // end loop over 't'
 
3796
        } // end loop over 's'
 
3797
      } // end loop over 'r'
3739
3798
      
3740
3799
      // Transform derivatives back to physical element
3741
3800
      for (unsigned int r = 0; r < num_derivatives; r++)
3743
3802
        for (unsigned int s = 0; s < num_derivatives; s++)
3744
3803
        {
3745
3804
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
3746
 
        }// end loop over 's'
3747
 
      }// end loop over 'r'
 
3805
        } // end loop over 's'
 
3806
      } // end loop over 'r'
3748
3807
        break;
3749
3808
      }
3750
3809
    case 10:
3799
3858
      for (unsigned int r = 0; r < 4; r++)
3800
3859
      {
3801
3860
        derivatives[r] = 0.0;
3802
 
      }// end loop over 'r'
 
3861
      } // end loop over 'r'
3803
3862
      
3804
3863
      // Declare derivative matrix (of polynomial basis).
3805
3864
      double dmats[6][6] = \
3833
3892
            dmats[t][u] = 1.0;
3834
3893
            }
3835
3894
            
3836
 
          }// end loop over 'u'
3837
 
        }// end loop over 't'
 
3895
          } // end loop over 'u'
 
3896
        } // end loop over 't'
3838
3897
        
3839
3898
        // Looping derivative order to generate dmats.
3840
3899
        for (unsigned int s = 0; s < n; s++)
3846
3905
            {
3847
3906
              dmats_old[t][u] = dmats[t][u];
3848
3907
              dmats[t][u] = 0.0;
3849
 
            }// end loop over 'u'
3850
 
          }// end loop over 't'
 
3908
            } // end loop over 'u'
 
3909
          } // end loop over 't'
3851
3910
          
3852
3911
          // Update dmats using an inner product.
3853
3912
          if (combinations[r][s] == 0)
3859
3918
              for (unsigned int tu = 0; tu < 6; tu++)
3860
3919
              {
3861
3920
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3862
 
              }// end loop over 'tu'
3863
 
            }// end loop over 'u'
3864
 
          }// end loop over 't'
 
3921
              } // end loop over 'tu'
 
3922
            } // end loop over 'u'
 
3923
          } // end loop over 't'
3865
3924
          }
3866
3925
          
3867
3926
          if (combinations[r][s] == 1)
3873
3932
              for (unsigned int tu = 0; tu < 6; tu++)
3874
3933
              {
3875
3934
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3876
 
              }// end loop over 'tu'
3877
 
            }// end loop over 'u'
3878
 
          }// end loop over 't'
 
3935
              } // end loop over 'tu'
 
3936
            } // end loop over 'u'
 
3937
          } // end loop over 't'
3879
3938
          }
3880
3939
          
3881
 
        }// end loop over 's'
 
3940
        } // end loop over 's'
3882
3941
        for (unsigned int s = 0; s < 6; s++)
3883
3942
        {
3884
3943
          for (unsigned int t = 0; t < 6; t++)
3885
3944
          {
3886
3945
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3887
 
          }// end loop over 't'
3888
 
        }// end loop over 's'
3889
 
      }// end loop over 'r'
 
3946
          } // end loop over 't'
 
3947
        } // end loop over 's'
 
3948
      } // end loop over 'r'
3890
3949
      
3891
3950
      // Transform derivatives back to physical element
3892
3951
      for (unsigned int r = 0; r < num_derivatives; r++)
3894
3953
        for (unsigned int s = 0; s < num_derivatives; s++)
3895
3954
        {
3896
3955
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
3897
 
        }// end loop over 's'
3898
 
      }// end loop over 'r'
 
3956
        } // end loop over 's'
 
3957
      } // end loop over 'r'
3899
3958
        break;
3900
3959
      }
3901
3960
    case 11:
3950
4009
      for (unsigned int r = 0; r < 4; r++)
3951
4010
      {
3952
4011
        derivatives[r] = 0.0;
3953
 
      }// end loop over 'r'
 
4012
      } // end loop over 'r'
3954
4013
      
3955
4014
      // Declare derivative matrix (of polynomial basis).
3956
4015
      double dmats[6][6] = \
3984
4043
            dmats[t][u] = 1.0;
3985
4044
            }
3986
4045
            
3987
 
          }// end loop over 'u'
3988
 
        }// end loop over 't'
 
4046
          } // end loop over 'u'
 
4047
        } // end loop over 't'
3989
4048
        
3990
4049
        // Looping derivative order to generate dmats.
3991
4050
        for (unsigned int s = 0; s < n; s++)
3997
4056
            {
3998
4057
              dmats_old[t][u] = dmats[t][u];
3999
4058
              dmats[t][u] = 0.0;
4000
 
            }// end loop over 'u'
4001
 
          }// end loop over 't'
 
4059
            } // end loop over 'u'
 
4060
          } // end loop over 't'
4002
4061
          
4003
4062
          // Update dmats using an inner product.
4004
4063
          if (combinations[r][s] == 0)
4010
4069
              for (unsigned int tu = 0; tu < 6; tu++)
4011
4070
              {
4012
4071
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4013
 
              }// end loop over 'tu'
4014
 
            }// end loop over 'u'
4015
 
          }// end loop over 't'
 
4072
              } // end loop over 'tu'
 
4073
            } // end loop over 'u'
 
4074
          } // end loop over 't'
4016
4075
          }
4017
4076
          
4018
4077
          if (combinations[r][s] == 1)
4024
4083
              for (unsigned int tu = 0; tu < 6; tu++)
4025
4084
              {
4026
4085
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4027
 
              }// end loop over 'tu'
4028
 
            }// end loop over 'u'
4029
 
          }// end loop over 't'
 
4086
              } // end loop over 'tu'
 
4087
            } // end loop over 'u'
 
4088
          } // end loop over 't'
4030
4089
          }
4031
4090
          
4032
 
        }// end loop over 's'
 
4091
        } // end loop over 's'
4033
4092
        for (unsigned int s = 0; s < 6; s++)
4034
4093
        {
4035
4094
          for (unsigned int t = 0; t < 6; t++)
4036
4095
          {
4037
4096
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4038
 
          }// end loop over 't'
4039
 
        }// end loop over 's'
4040
 
      }// end loop over 'r'
 
4097
          } // end loop over 't'
 
4098
        } // end loop over 's'
 
4099
      } // end loop over 'r'
4041
4100
      
4042
4101
      // Transform derivatives back to physical element
4043
4102
      for (unsigned int r = 0; r < num_derivatives; r++)
4045
4104
        for (unsigned int s = 0; s < num_derivatives; s++)
4046
4105
        {
4047
4106
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
4048
 
        }// end loop over 's'
4049
 
      }// end loop over 'r'
 
4107
        } // end loop over 's'
 
4108
      } // end loop over 'r'
4050
4109
        break;
4051
4110
      }
4052
4111
    }
4053
4112
    
4054
4113
  }
4055
4114
 
4056
 
  /// Evaluate order n derivatives of all basis functions at given point x in cell
4057
 
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
4115
  /// Evaluate order n derivatives of basis function i at given point x in cell (non-static member function)
 
4116
  virtual void evaluate_basis_derivatives(std::size_t i,
 
4117
                                          std::size_t n,
 
4118
                                          double* values,
 
4119
                                          const double* x,
 
4120
                                          const double* vertex_coordinates,
 
4121
                                          int cell_orientation) const
 
4122
  {
 
4123
    _evaluate_basis_derivatives(i, n, values, x, vertex_coordinates, cell_orientation);
 
4124
  }
 
4125
 
 
4126
  /// Evaluate order n derivatives of all basis functions at given point x in cell (actual implementation)
 
4127
  static void _evaluate_basis_derivatives_all(std::size_t n,
4058
4128
                                              double* values,
4059
4129
                                              const double* x,
4060
4130
                                              const double* vertex_coordinates,
4061
 
                                              int cell_orientation) const
 
4131
                                              int cell_orientation)
4062
4132
  {
4063
4133
    // Call evaluate_basis_all if order of derivatives is equal to zero.
4064
4134
    if (n == 0)
4065
4135
    {
4066
 
      evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
4136
      _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
4067
4137
      return ;
4068
4138
    }
4069
4139
    
4072
4142
    for (unsigned int r = 0; r < n; r++)
4073
4143
    {
4074
4144
      num_derivatives *= 2;
4075
 
    }// end loop over 'r'
 
4145
    } // end loop over 'r'
4076
4146
    
4077
4147
    // Set values equal to zero.
4078
4148
    for (unsigned int r = 0; r < 12; r++)
4080
4150
      for (unsigned int s = 0; s < 2*num_derivatives; s++)
4081
4151
      {
4082
4152
        values[r*2*num_derivatives + s] = 0.0;
4083
 
      }// end loop over 's'
4084
 
    }// end loop over 'r'
 
4153
      } // end loop over 's'
 
4154
    } // end loop over 'r'
4085
4155
    
4086
4156
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
4087
4157
    if (n > 2)
4094
4164
    for (unsigned int r = 0; r < 8; r++)
4095
4165
    {
4096
4166
      dof_values[r] = 0.0;
4097
 
    }// end loop over 'r'
 
4167
    } // end loop over 'r'
4098
4168
    
4099
4169
    // Loop dofs and call evaluate_basis_derivatives.
4100
4170
    for (unsigned int r = 0; r < 12; r++)
4101
4171
    {
4102
 
      evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
 
4172
      _evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
4103
4173
      for (unsigned int s = 0; s < 2*num_derivatives; s++)
4104
4174
      {
4105
4175
        values[r*2*num_derivatives + s] = dof_values[s];
4106
 
      }// end loop over 's'
4107
 
    }// end loop over 'r'
 
4176
      } // end loop over 's'
 
4177
    } // end loop over 'r'
 
4178
  }
 
4179
 
 
4180
  /// Evaluate order n derivatives of all basis functions at given point x in cell (non-static member function)
 
4181
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
4182
                                              double* values,
 
4183
                                              const double* x,
 
4184
                                              const double* vertex_coordinates,
 
4185
                                              int cell_orientation) const
 
4186
  {
 
4187
    _evaluate_basis_derivatives_all(n, values, x, vertex_coordinates, cell_orientation);
4108
4188
  }
4109
4189
 
4110
4190
  /// Evaluate linear functional for dof i on the function f
4372
4452
  /// Return a string identifying the finite element
4373
4453
  virtual const char* signature() const
4374
4454
  {
4375
 
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
 
4455
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 1, None)";
4376
4456
  }
4377
4457
 
4378
4458
  /// Return the cell shape
4411
4491
    return 1;
4412
4492
  }
4413
4493
 
4414
 
  /// Evaluate basis function i at given point x in cell
4415
 
  virtual void evaluate_basis(std::size_t i,
 
4494
  /// Evaluate basis function i at given point x in cell (actual implementation)
 
4495
  static void _evaluate_basis(std::size_t i,
4416
4496
                              double* values,
4417
4497
                              const double* x,
4418
4498
                              const double* vertex_coordinates,
4419
 
                              int cell_orientation) const
 
4499
                              int cell_orientation)
4420
4500
  {
4421
4501
    // Compute Jacobian
4422
4502
    double J[4];
4465
4545
      for (unsigned int r = 0; r < 3; r++)
4466
4546
      {
4467
4547
        *values += coefficients0[r]*basisvalues[r];
4468
 
      }// end loop over 'r'
 
4548
      } // end loop over 'r'
4469
4549
        break;
4470
4550
      }
4471
4551
    case 1:
4493
4573
      for (unsigned int r = 0; r < 3; r++)
4494
4574
      {
4495
4575
        *values += coefficients0[r]*basisvalues[r];
4496
 
      }// end loop over 'r'
 
4576
      } // end loop over 'r'
4497
4577
        break;
4498
4578
      }
4499
4579
    case 2:
4521
4601
      for (unsigned int r = 0; r < 3; r++)
4522
4602
      {
4523
4603
        *values += coefficients0[r]*basisvalues[r];
4524
 
      }// end loop over 'r'
 
4604
      } // end loop over 'r'
4525
4605
        break;
4526
4606
      }
4527
4607
    }
4528
4608
    
4529
4609
  }
4530
4610
 
4531
 
  /// Evaluate all basis functions at given point x in cell
4532
 
  virtual void evaluate_basis_all(double* values,
 
4611
  /// Evaluate basis function i at given point x in cell (non-static member function)
 
4612
  virtual void evaluate_basis(std::size_t i,
 
4613
                              double* values,
 
4614
                              const double* x,
 
4615
                              const double* vertex_coordinates,
 
4616
                              int cell_orientation) const
 
4617
  {
 
4618
    _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
 
4619
  }
 
4620
 
 
4621
  /// Evaluate all basis functions at given point x in cell (actual implementation)
 
4622
  static void _evaluate_basis_all(double* values,
4533
4623
                                  const double* x,
4534
4624
                                  const double* vertex_coordinates,
4535
 
                                  int cell_orientation) const
 
4625
                                  int cell_orientation)
4536
4626
  {
4537
4627
    // Helper variable to hold values of a single dof.
4538
4628
    double dof_values = 0.0;
4540
4630
    // Loop dofs and call evaluate_basis
4541
4631
    for (unsigned int r = 0; r < 3; r++)
4542
4632
    {
4543
 
      evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
 
4633
      _evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
4544
4634
      values[r] = dof_values;
4545
 
    }// end loop over 'r'
4546
 
  }
4547
 
 
4548
 
  /// Evaluate order n derivatives of basis function i at given point x in cell
4549
 
  virtual void evaluate_basis_derivatives(std::size_t i,
 
4635
    } // end loop over 'r'
 
4636
  }
 
4637
 
 
4638
  /// Evaluate all basis functions at given point x in cell (non-static member function)
 
4639
  virtual void evaluate_basis_all(double* values,
 
4640
                                  const double* x,
 
4641
                                  const double* vertex_coordinates,
 
4642
                                  int cell_orientation) const
 
4643
  {
 
4644
    _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
4645
  }
 
4646
 
 
4647
  /// Evaluate order n derivatives of basis function i at given point x in cell (actual implementation)
 
4648
  static void _evaluate_basis_derivatives(std::size_t i,
4550
4649
                                          std::size_t n,
4551
4650
                                          double* values,
4552
4651
                                          const double* x,
4553
4652
                                          const double* vertex_coordinates,
4554
 
                                          int cell_orientation) const
 
4653
                                          int cell_orientation)
4555
4654
  {
4556
4655
    
4557
4656
    // Compute number of derivatives.
4559
4658
    for (unsigned int r = 0; r < n; r++)
4560
4659
    {
4561
4660
      num_derivatives *= 2;
4562
 
    }// end loop over 'r'
 
4661
    } // end loop over 'r'
4563
4662
    
4564
4663
    // Reset values. Assuming that values is always an array.
4565
4664
    for (unsigned int r = 0; r < num_derivatives; r++)
4566
4665
    {
4567
4666
      values[r] = 0.0;
4568
 
    }// end loop over 'r'
 
4667
    } // end loop over 'r'
4569
4668
    
4570
4669
    // Call evaluate_basis if order of derivatives is equal to zero.
4571
4670
    if (n == 0)
4572
4671
    {
4573
 
      evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
 
4672
      _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
4574
4673
      return ;
4575
4674
    }
4576
4675
    
4685
4784
      for (unsigned int r = 0; r < 2; r++)
4686
4785
      {
4687
4786
        derivatives[r] = 0.0;
4688
 
      }// end loop over 'r'
 
4787
      } // end loop over 'r'
4689
4788
      
4690
4789
      // Declare derivative matrix (of polynomial basis).
4691
4790
      double dmats[3][3] = \
4713
4812
            dmats[t][u] = 1.0;
4714
4813
            }
4715
4814
            
4716
 
          }// end loop over 'u'
4717
 
        }// end loop over 't'
 
4815
          } // end loop over 'u'
 
4816
        } // end loop over 't'
4718
4817
        
4719
4818
        // Looping derivative order to generate dmats.
4720
4819
        for (unsigned int s = 0; s < n; s++)
4726
4825
            {
4727
4826
              dmats_old[t][u] = dmats[t][u];
4728
4827
              dmats[t][u] = 0.0;
4729
 
            }// end loop over 'u'
4730
 
          }// end loop over 't'
 
4828
            } // end loop over 'u'
 
4829
          } // end loop over 't'
4731
4830
          
4732
4831
          // Update dmats using an inner product.
4733
4832
          if (combinations[r][s] == 0)
4739
4838
              for (unsigned int tu = 0; tu < 3; tu++)
4740
4839
              {
4741
4840
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4742
 
              }// end loop over 'tu'
4743
 
            }// end loop over 'u'
4744
 
          }// end loop over 't'
 
4841
              } // end loop over 'tu'
 
4842
            } // end loop over 'u'
 
4843
          } // end loop over 't'
4745
4844
          }
4746
4845
          
4747
4846
          if (combinations[r][s] == 1)
4753
4852
              for (unsigned int tu = 0; tu < 3; tu++)
4754
4853
              {
4755
4854
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4756
 
              }// end loop over 'tu'
4757
 
            }// end loop over 'u'
4758
 
          }// end loop over 't'
 
4855
              } // end loop over 'tu'
 
4856
            } // end loop over 'u'
 
4857
          } // end loop over 't'
4759
4858
          }
4760
4859
          
4761
 
        }// end loop over 's'
 
4860
        } // end loop over 's'
4762
4861
        for (unsigned int s = 0; s < 3; s++)
4763
4862
        {
4764
4863
          for (unsigned int t = 0; t < 3; t++)
4765
4864
          {
4766
4865
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4767
 
          }// end loop over 't'
4768
 
        }// end loop over 's'
4769
 
      }// end loop over 'r'
 
4866
          } // end loop over 't'
 
4867
        } // end loop over 's'
 
4868
      } // end loop over 'r'
4770
4869
      
4771
4870
      // Transform derivatives back to physical element
4772
4871
      for (unsigned int r = 0; r < num_derivatives; r++)
4774
4873
        for (unsigned int s = 0; s < num_derivatives; s++)
4775
4874
        {
4776
4875
          values[r] += transform[r][s]*derivatives[s];
4777
 
        }// end loop over 's'
4778
 
      }// end loop over 'r'
 
4876
        } // end loop over 's'
 
4877
      } // end loop over 'r'
4779
4878
        break;
4780
4879
      }
4781
4880
    case 1:
4816
4915
      for (unsigned int r = 0; r < 2; r++)
4817
4916
      {
4818
4917
        derivatives[r] = 0.0;
4819
 
      }// end loop over 'r'
 
4918
      } // end loop over 'r'
4820
4919
      
4821
4920
      // Declare derivative matrix (of polynomial basis).
4822
4921
      double dmats[3][3] = \
4844
4943
            dmats[t][u] = 1.0;
4845
4944
            }
4846
4945
            
4847
 
          }// end loop over 'u'
4848
 
        }// end loop over 't'
 
4946
          } // end loop over 'u'
 
4947
        } // end loop over 't'
4849
4948
        
4850
4949
        // Looping derivative order to generate dmats.
4851
4950
        for (unsigned int s = 0; s < n; s++)
4857
4956
            {
4858
4957
              dmats_old[t][u] = dmats[t][u];
4859
4958
              dmats[t][u] = 0.0;
4860
 
            }// end loop over 'u'
4861
 
          }// end loop over 't'
 
4959
            } // end loop over 'u'
 
4960
          } // end loop over 't'
4862
4961
          
4863
4962
          // Update dmats using an inner product.
4864
4963
          if (combinations[r][s] == 0)
4870
4969
              for (unsigned int tu = 0; tu < 3; tu++)
4871
4970
              {
4872
4971
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4873
 
              }// end loop over 'tu'
4874
 
            }// end loop over 'u'
4875
 
          }// end loop over 't'
 
4972
              } // end loop over 'tu'
 
4973
            } // end loop over 'u'
 
4974
          } // end loop over 't'
4876
4975
          }
4877
4976
          
4878
4977
          if (combinations[r][s] == 1)
4884
4983
              for (unsigned int tu = 0; tu < 3; tu++)
4885
4984
              {
4886
4985
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4887
 
              }// end loop over 'tu'
4888
 
            }// end loop over 'u'
4889
 
          }// end loop over 't'
 
4986
              } // end loop over 'tu'
 
4987
            } // end loop over 'u'
 
4988
          } // end loop over 't'
4890
4989
          }
4891
4990
          
4892
 
        }// end loop over 's'
 
4991
        } // end loop over 's'
4893
4992
        for (unsigned int s = 0; s < 3; s++)
4894
4993
        {
4895
4994
          for (unsigned int t = 0; t < 3; t++)
4896
4995
          {
4897
4996
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4898
 
          }// end loop over 't'
4899
 
        }// end loop over 's'
4900
 
      }// end loop over 'r'
 
4997
          } // end loop over 't'
 
4998
        } // end loop over 's'
 
4999
      } // end loop over 'r'
4901
5000
      
4902
5001
      // Transform derivatives back to physical element
4903
5002
      for (unsigned int r = 0; r < num_derivatives; r++)
4905
5004
        for (unsigned int s = 0; s < num_derivatives; s++)
4906
5005
        {
4907
5006
          values[r] += transform[r][s]*derivatives[s];
4908
 
        }// end loop over 's'
4909
 
      }// end loop over 'r'
 
5007
        } // end loop over 's'
 
5008
      } // end loop over 'r'
4910
5009
        break;
4911
5010
      }
4912
5011
    case 2:
4947
5046
      for (unsigned int r = 0; r < 2; r++)
4948
5047
      {
4949
5048
        derivatives[r] = 0.0;
4950
 
      }// end loop over 'r'
 
5049
      } // end loop over 'r'
4951
5050
      
4952
5051
      // Declare derivative matrix (of polynomial basis).
4953
5052
      double dmats[3][3] = \
4975
5074
            dmats[t][u] = 1.0;
4976
5075
            }
4977
5076
            
4978
 
          }// end loop over 'u'
4979
 
        }// end loop over 't'
 
5077
          } // end loop over 'u'
 
5078
        } // end loop over 't'
4980
5079
        
4981
5080
        // Looping derivative order to generate dmats.
4982
5081
        for (unsigned int s = 0; s < n; s++)
4988
5087
            {
4989
5088
              dmats_old[t][u] = dmats[t][u];
4990
5089
              dmats[t][u] = 0.0;
4991
 
            }// end loop over 'u'
4992
 
          }// end loop over 't'
 
5090
            } // end loop over 'u'
 
5091
          } // end loop over 't'
4993
5092
          
4994
5093
          // Update dmats using an inner product.
4995
5094
          if (combinations[r][s] == 0)
5001
5100
              for (unsigned int tu = 0; tu < 3; tu++)
5002
5101
              {
5003
5102
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
5004
 
              }// end loop over 'tu'
5005
 
            }// end loop over 'u'
5006
 
          }// end loop over 't'
 
5103
              } // end loop over 'tu'
 
5104
            } // end loop over 'u'
 
5105
          } // end loop over 't'
5007
5106
          }
5008
5107
          
5009
5108
          if (combinations[r][s] == 1)
5015
5114
              for (unsigned int tu = 0; tu < 3; tu++)
5016
5115
              {
5017
5116
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
5018
 
              }// end loop over 'tu'
5019
 
            }// end loop over 'u'
5020
 
          }// end loop over 't'
 
5117
              } // end loop over 'tu'
 
5118
            } // end loop over 'u'
 
5119
          } // end loop over 't'
5021
5120
          }
5022
5121
          
5023
 
        }// end loop over 's'
 
5122
        } // end loop over 's'
5024
5123
        for (unsigned int s = 0; s < 3; s++)
5025
5124
        {
5026
5125
          for (unsigned int t = 0; t < 3; t++)
5027
5126
          {
5028
5127
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
5029
 
          }// end loop over 't'
5030
 
        }// end loop over 's'
5031
 
      }// end loop over 'r'
 
5128
          } // end loop over 't'
 
5129
        } // end loop over 's'
 
5130
      } // end loop over 'r'
5032
5131
      
5033
5132
      // Transform derivatives back to physical element
5034
5133
      for (unsigned int r = 0; r < num_derivatives; r++)
5036
5135
        for (unsigned int s = 0; s < num_derivatives; s++)
5037
5136
        {
5038
5137
          values[r] += transform[r][s]*derivatives[s];
5039
 
        }// end loop over 's'
5040
 
      }// end loop over 'r'
 
5138
        } // end loop over 's'
 
5139
      } // end loop over 'r'
5041
5140
        break;
5042
5141
      }
5043
5142
    }
5044
5143
    
5045
5144
  }
5046
5145
 
5047
 
  /// Evaluate order n derivatives of all basis functions at given point x in cell
5048
 
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
5146
  /// Evaluate order n derivatives of basis function i at given point x in cell (non-static member function)
 
5147
  virtual void evaluate_basis_derivatives(std::size_t i,
 
5148
                                          std::size_t n,
 
5149
                                          double* values,
 
5150
                                          const double* x,
 
5151
                                          const double* vertex_coordinates,
 
5152
                                          int cell_orientation) const
 
5153
  {
 
5154
    _evaluate_basis_derivatives(i, n, values, x, vertex_coordinates, cell_orientation);
 
5155
  }
 
5156
 
 
5157
  /// Evaluate order n derivatives of all basis functions at given point x in cell (actual implementation)
 
5158
  static void _evaluate_basis_derivatives_all(std::size_t n,
5049
5159
                                              double* values,
5050
5160
                                              const double* x,
5051
5161
                                              const double* vertex_coordinates,
5052
 
                                              int cell_orientation) const
 
5162
                                              int cell_orientation)
5053
5163
  {
5054
5164
    // Call evaluate_basis_all if order of derivatives is equal to zero.
5055
5165
    if (n == 0)
5056
5166
    {
5057
 
      evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
5167
      _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
5058
5168
      return ;
5059
5169
    }
5060
5170
    
5063
5173
    for (unsigned int r = 0; r < n; r++)
5064
5174
    {
5065
5175
      num_derivatives *= 2;
5066
 
    }// end loop over 'r'
 
5176
    } // end loop over 'r'
5067
5177
    
5068
5178
    // Set values equal to zero.
5069
5179
    for (unsigned int r = 0; r < 3; r++)
5071
5181
      for (unsigned int s = 0; s < num_derivatives; s++)
5072
5182
      {
5073
5183
        values[r*num_derivatives + s] = 0.0;
5074
 
      }// end loop over 's'
5075
 
    }// end loop over 'r'
 
5184
      } // end loop over 's'
 
5185
    } // end loop over 'r'
5076
5186
    
5077
5187
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
5078
5188
    if (n > 1)
5085
5195
    for (unsigned int r = 0; r < 2; r++)
5086
5196
    {
5087
5197
      dof_values[r] = 0.0;
5088
 
    }// end loop over 'r'
 
5198
    } // end loop over 'r'
5089
5199
    
5090
5200
    // Loop dofs and call evaluate_basis_derivatives.
5091
5201
    for (unsigned int r = 0; r < 3; r++)
5092
5202
    {
5093
 
      evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
 
5203
      _evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
5094
5204
      for (unsigned int s = 0; s < num_derivatives; s++)
5095
5205
      {
5096
5206
        values[r*num_derivatives + s] = dof_values[s];
5097
 
      }// end loop over 's'
5098
 
    }// end loop over 'r'
 
5207
      } // end loop over 's'
 
5208
    } // end loop over 'r'
 
5209
  }
 
5210
 
 
5211
  /// Evaluate order n derivatives of all basis functions at given point x in cell (non-static member function)
 
5212
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
5213
                                              double* values,
 
5214
                                              const double* x,
 
5215
                                              const double* vertex_coordinates,
 
5216
                                              int cell_orientation) const
 
5217
  {
 
5218
    _evaluate_basis_derivatives_all(n, values, x, vertex_coordinates, cell_orientation);
5099
5219
  }
5100
5220
 
5101
5221
  /// Evaluate linear functional for dof i on the function f
5237
5357
  /// Return a string identifying the finite element
5238
5358
  virtual const char* signature() const
5239
5359
  {
5240
 
    return "MixedElement(*[VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (3,) })";
 
5360
    return "MixedElement(VectorElement('Lagrange', Domain(Cell('triangle', 2)), 2, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 1, None), **{'value_shape': (3,) })";
5241
5361
  }
5242
5362
 
5243
5363
  /// Return the cell shape
5285
5405
    return 0;
5286
5406
  }
5287
5407
 
5288
 
  /// Evaluate basis function i at given point x in cell
5289
 
  virtual void evaluate_basis(std::size_t i,
 
5408
  /// Evaluate basis function i at given point x in cell (actual implementation)
 
5409
  static void _evaluate_basis(std::size_t i,
5290
5410
                              double* values,
5291
5411
                              const double* x,
5292
5412
                              const double* vertex_coordinates,
5293
 
                              int cell_orientation) const
 
5413
                              int cell_orientation)
5294
5414
  {
5295
5415
    // Compute Jacobian
5296
5416
    double J[4];
5349
5469
      for (unsigned int r = 0; r < 6; r++)
5350
5470
      {
5351
5471
        values[0] += coefficients0[r]*basisvalues[r];
5352
 
      }// end loop over 'r'
 
5472
      } // end loop over 'r'
5353
5473
        break;
5354
5474
      }
5355
5475
    case 1:
5385
5505
      for (unsigned int r = 0; r < 6; r++)
5386
5506
      {
5387
5507
        values[0] += coefficients0[r]*basisvalues[r];
5388
 
      }// end loop over 'r'
 
5508
      } // end loop over 'r'
5389
5509
        break;
5390
5510
      }
5391
5511
    case 2:
5421
5541
      for (unsigned int r = 0; r < 6; r++)
5422
5542
      {
5423
5543
        values[0] += coefficients0[r]*basisvalues[r];
5424
 
      }// end loop over 'r'
 
5544
      } // end loop over 'r'
5425
5545
        break;
5426
5546
      }
5427
5547
    case 3:
5457
5577
      for (unsigned int r = 0; r < 6; r++)
5458
5578
      {
5459
5579
        values[0] += coefficients0[r]*basisvalues[r];
5460
 
      }// end loop over 'r'
 
5580
      } // end loop over 'r'
5461
5581
        break;
5462
5582
      }
5463
5583
    case 4:
5493
5613
      for (unsigned int r = 0; r < 6; r++)
5494
5614
      {
5495
5615
        values[0] += coefficients0[r]*basisvalues[r];
5496
 
      }// end loop over 'r'
 
5616
      } // end loop over 'r'
5497
5617
        break;
5498
5618
      }
5499
5619
    case 5:
5529
5649
      for (unsigned int r = 0; r < 6; r++)
5530
5650
      {
5531
5651
        values[0] += coefficients0[r]*basisvalues[r];
5532
 
      }// end loop over 'r'
 
5652
      } // end loop over 'r'
5533
5653
        break;
5534
5654
      }
5535
5655
    case 6:
5565
5685
      for (unsigned int r = 0; r < 6; r++)
5566
5686
      {
5567
5687
        values[1] += coefficients0[r]*basisvalues[r];
5568
 
      }// end loop over 'r'
 
5688
      } // end loop over 'r'
5569
5689
        break;
5570
5690
      }
5571
5691
    case 7:
5601
5721
      for (unsigned int r = 0; r < 6; r++)
5602
5722
      {
5603
5723
        values[1] += coefficients0[r]*basisvalues[r];
5604
 
      }// end loop over 'r'
 
5724
      } // end loop over 'r'
5605
5725
        break;
5606
5726
      }
5607
5727
    case 8:
5637
5757
      for (unsigned int r = 0; r < 6; r++)
5638
5758
      {
5639
5759
        values[1] += coefficients0[r]*basisvalues[r];
5640
 
      }// end loop over 'r'
 
5760
      } // end loop over 'r'
5641
5761
        break;
5642
5762
      }
5643
5763
    case 9:
5673
5793
      for (unsigned int r = 0; r < 6; r++)
5674
5794
      {
5675
5795
        values[1] += coefficients0[r]*basisvalues[r];
5676
 
      }// end loop over 'r'
 
5796
      } // end loop over 'r'
5677
5797
        break;
5678
5798
      }
5679
5799
    case 10:
5709
5829
      for (unsigned int r = 0; r < 6; r++)
5710
5830
      {
5711
5831
        values[1] += coefficients0[r]*basisvalues[r];
5712
 
      }// end loop over 'r'
 
5832
      } // end loop over 'r'
5713
5833
        break;
5714
5834
      }
5715
5835
    case 11:
5745
5865
      for (unsigned int r = 0; r < 6; r++)
5746
5866
      {
5747
5867
        values[1] += coefficients0[r]*basisvalues[r];
5748
 
      }// end loop over 'r'
 
5868
      } // end loop over 'r'
5749
5869
        break;
5750
5870
      }
5751
5871
    case 12:
5773
5893
      for (unsigned int r = 0; r < 3; r++)
5774
5894
      {
5775
5895
        values[2] += coefficients0[r]*basisvalues[r];
5776
 
      }// end loop over 'r'
 
5896
      } // end loop over 'r'
5777
5897
        break;
5778
5898
      }
5779
5899
    case 13:
5801
5921
      for (unsigned int r = 0; r < 3; r++)
5802
5922
      {
5803
5923
        values[2] += coefficients0[r]*basisvalues[r];
5804
 
      }// end loop over 'r'
 
5924
      } // end loop over 'r'
5805
5925
        break;
5806
5926
      }
5807
5927
    case 14:
5829
5949
      for (unsigned int r = 0; r < 3; r++)
5830
5950
      {
5831
5951
        values[2] += coefficients0[r]*basisvalues[r];
5832
 
      }// end loop over 'r'
 
5952
      } // end loop over 'r'
5833
5953
        break;
5834
5954
      }
5835
5955
    }
5836
5956
    
5837
5957
  }
5838
5958
 
5839
 
  /// Evaluate all basis functions at given point x in cell
5840
 
  virtual void evaluate_basis_all(double* values,
 
5959
  /// Evaluate basis function i at given point x in cell (non-static member function)
 
5960
  virtual void evaluate_basis(std::size_t i,
 
5961
                              double* values,
 
5962
                              const double* x,
 
5963
                              const double* vertex_coordinates,
 
5964
                              int cell_orientation) const
 
5965
  {
 
5966
    _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
 
5967
  }
 
5968
 
 
5969
  /// Evaluate all basis functions at given point x in cell (actual implementation)
 
5970
  static void _evaluate_basis_all(double* values,
5841
5971
                                  const double* x,
5842
5972
                                  const double* vertex_coordinates,
5843
 
                                  int cell_orientation) const
 
5973
                                  int cell_orientation)
5844
5974
  {
5845
5975
    // Helper variable to hold values of a single dof.
5846
5976
    double dof_values[3] = {0.0, 0.0, 0.0};
5848
5978
    // Loop dofs and call evaluate_basis
5849
5979
    for (unsigned int r = 0; r < 15; r++)
5850
5980
    {
5851
 
      evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
 
5981
      _evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
5852
5982
      for (unsigned int s = 0; s < 3; s++)
5853
5983
      {
5854
5984
        values[r*3 + s] = dof_values[s];
5855
 
      }// end loop over 's'
5856
 
    }// end loop over 'r'
5857
 
  }
5858
 
 
5859
 
  /// Evaluate order n derivatives of basis function i at given point x in cell
5860
 
  virtual void evaluate_basis_derivatives(std::size_t i,
 
5985
      } // end loop over 's'
 
5986
    } // end loop over 'r'
 
5987
  }
 
5988
 
 
5989
  /// Evaluate all basis functions at given point x in cell (non-static member function)
 
5990
  virtual void evaluate_basis_all(double* values,
 
5991
                                  const double* x,
 
5992
                                  const double* vertex_coordinates,
 
5993
                                  int cell_orientation) const
 
5994
  {
 
5995
    _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
5996
  }
 
5997
 
 
5998
  /// Evaluate order n derivatives of basis function i at given point x in cell (actual implementation)
 
5999
  static void _evaluate_basis_derivatives(std::size_t i,
5861
6000
                                          std::size_t n,
5862
6001
                                          double* values,
5863
6002
                                          const double* x,
5864
6003
                                          const double* vertex_coordinates,
5865
 
                                          int cell_orientation) const
 
6004
                                          int cell_orientation)
5866
6005
  {
5867
6006
    
5868
6007
    // Compute number of derivatives.
5870
6009
    for (unsigned int r = 0; r < n; r++)
5871
6010
    {
5872
6011
      num_derivatives *= 2;
5873
 
    }// end loop over 'r'
 
6012
    } // end loop over 'r'
5874
6013
    
5875
6014
    // Reset values. Assuming that values is always an array.
5876
6015
    for (unsigned int r = 0; r < 3*num_derivatives; r++)
5877
6016
    {
5878
6017
      values[r] = 0.0;
5879
 
    }// end loop over 'r'
 
6018
    } // end loop over 'r'
5880
6019
    
5881
6020
    // Call evaluate_basis if order of derivatives is equal to zero.
5882
6021
    if (n == 0)
5883
6022
    {
5884
 
      evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
 
6023
      _evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
5885
6024
      return ;
5886
6025
    }
5887
6026
    
6010
6149
      for (unsigned int r = 0; r < 4; r++)
6011
6150
      {
6012
6151
        derivatives[r] = 0.0;
6013
 
      }// end loop over 'r'
 
6152
      } // end loop over 'r'
6014
6153
      
6015
6154
      // Declare derivative matrix (of polynomial basis).
6016
6155
      double dmats[6][6] = \
6044
6183
            dmats[t][u] = 1.0;
6045
6184
            }
6046
6185
            
6047
 
          }// end loop over 'u'
6048
 
        }// end loop over 't'
 
6186
          } // end loop over 'u'
 
6187
        } // end loop over 't'
6049
6188
        
6050
6189
        // Looping derivative order to generate dmats.
6051
6190
        for (unsigned int s = 0; s < n; s++)
6057
6196
            {
6058
6197
              dmats_old[t][u] = dmats[t][u];
6059
6198
              dmats[t][u] = 0.0;
6060
 
            }// end loop over 'u'
6061
 
          }// end loop over 't'
 
6199
            } // end loop over 'u'
 
6200
          } // end loop over 't'
6062
6201
          
6063
6202
          // Update dmats using an inner product.
6064
6203
          if (combinations[r][s] == 0)
6070
6209
              for (unsigned int tu = 0; tu < 6; tu++)
6071
6210
              {
6072
6211
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
6073
 
              }// end loop over 'tu'
6074
 
            }// end loop over 'u'
6075
 
          }// end loop over 't'
 
6212
              } // end loop over 'tu'
 
6213
            } // end loop over 'u'
 
6214
          } // end loop over 't'
6076
6215
          }
6077
6216
          
6078
6217
          if (combinations[r][s] == 1)
6084
6223
              for (unsigned int tu = 0; tu < 6; tu++)
6085
6224
              {
6086
6225
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
6087
 
              }// end loop over 'tu'
6088
 
            }// end loop over 'u'
6089
 
          }// end loop over 't'
 
6226
              } // end loop over 'tu'
 
6227
            } // end loop over 'u'
 
6228
          } // end loop over 't'
6090
6229
          }
6091
6230
          
6092
 
        }// end loop over 's'
 
6231
        } // end loop over 's'
6093
6232
        for (unsigned int s = 0; s < 6; s++)
6094
6233
        {
6095
6234
          for (unsigned int t = 0; t < 6; t++)
6096
6235
          {
6097
6236
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
6098
 
          }// end loop over 't'
6099
 
        }// end loop over 's'
6100
 
      }// end loop over 'r'
 
6237
          } // end loop over 't'
 
6238
        } // end loop over 's'
 
6239
      } // end loop over 'r'
6101
6240
      
6102
6241
      // Transform derivatives back to physical element
6103
6242
      for (unsigned int r = 0; r < num_derivatives; r++)
6105
6244
        for (unsigned int s = 0; s < num_derivatives; s++)
6106
6245
        {
6107
6246
          values[r] += transform[r][s]*derivatives[s];
6108
 
        }// end loop over 's'
6109
 
      }// end loop over 'r'
 
6247
        } // end loop over 's'
 
6248
      } // end loop over 'r'
6110
6249
        break;
6111
6250
      }
6112
6251
    case 1:
6161
6300
      for (unsigned int r = 0; r < 4; r++)
6162
6301
      {
6163
6302
        derivatives[r] = 0.0;
6164
 
      }// end loop over 'r'
 
6303
      } // end loop over 'r'
6165
6304
      
6166
6305
      // Declare derivative matrix (of polynomial basis).
6167
6306
      double dmats[6][6] = \
6195
6334
            dmats[t][u] = 1.0;
6196
6335
            }
6197
6336
            
6198
 
          }// end loop over 'u'
6199
 
        }// end loop over 't'
 
6337
          } // end loop over 'u'
 
6338
        } // end loop over 't'
6200
6339
        
6201
6340
        // Looping derivative order to generate dmats.
6202
6341
        for (unsigned int s = 0; s < n; s++)
6208
6347
            {
6209
6348
              dmats_old[t][u] = dmats[t][u];
6210
6349
              dmats[t][u] = 0.0;
6211
 
            }// end loop over 'u'
6212
 
          }// end loop over 't'
 
6350
            } // end loop over 'u'
 
6351
          } // end loop over 't'
6213
6352
          
6214
6353
          // Update dmats using an inner product.
6215
6354
          if (combinations[r][s] == 0)
6221
6360
              for (unsigned int tu = 0; tu < 6; tu++)
6222
6361
              {
6223
6362
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
6224
 
              }// end loop over 'tu'
6225
 
            }// end loop over 'u'
6226
 
          }// end loop over 't'
 
6363
              } // end loop over 'tu'
 
6364
            } // end loop over 'u'
 
6365
          } // end loop over 't'
6227
6366
          }
6228
6367
          
6229
6368
          if (combinations[r][s] == 1)
6235
6374
              for (unsigned int tu = 0; tu < 6; tu++)
6236
6375
              {
6237
6376
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
6238
 
              }// end loop over 'tu'
6239
 
            }// end loop over 'u'
6240
 
          }// end loop over 't'
 
6377
              } // end loop over 'tu'
 
6378
            } // end loop over 'u'
 
6379
          } // end loop over 't'
6241
6380
          }
6242
6381
          
6243
 
        }// end loop over 's'
 
6382
        } // end loop over 's'
6244
6383
        for (unsigned int s = 0; s < 6; s++)
6245
6384
        {
6246
6385
          for (unsigned int t = 0; t < 6; t++)
6247
6386
          {
6248
6387
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
6249
 
          }// end loop over 't'
6250
 
        }// end loop over 's'
6251
 
      }// end loop over 'r'
 
6388
          } // end loop over 't'
 
6389
        } // end loop over 's'
 
6390
      } // end loop over 'r'
6252
6391
      
6253
6392
      // Transform derivatives back to physical element
6254
6393
      for (unsigned int r = 0; r < num_derivatives; r++)
6256
6395
        for (unsigned int s = 0; s < num_derivatives; s++)
6257
6396
        {
6258
6397
          values[r] += transform[r][s]*derivatives[s];
6259
 
        }// end loop over 's'
6260
 
      }// end loop over 'r'
 
6398
        } // end loop over 's'
 
6399
      } // end loop over 'r'
6261
6400
        break;
6262
6401
      }
6263
6402
    case 2:
6312
6451
      for (unsigned int r = 0; r < 4; r++)
6313
6452
      {
6314
6453
        derivatives[r] = 0.0;
6315
 
      }// end loop over 'r'
 
6454
      } // end loop over 'r'
6316
6455
      
6317
6456
      // Declare derivative matrix (of polynomial basis).
6318
6457
      double dmats[6][6] = \
6346
6485
            dmats[t][u] = 1.0;
6347
6486
            }
6348
6487
            
6349
 
          }// end loop over 'u'
6350
 
        }// end loop over 't'
 
6488
          } // end loop over 'u'
 
6489
        } // end loop over 't'
6351
6490
        
6352
6491
        // Looping derivative order to generate dmats.
6353
6492
        for (unsigned int s = 0; s < n; s++)
6359
6498
            {
6360
6499
              dmats_old[t][u] = dmats[t][u];
6361
6500
              dmats[t][u] = 0.0;
6362
 
            }// end loop over 'u'
6363
 
          }// end loop over 't'
 
6501
            } // end loop over 'u'
 
6502
          } // end loop over 't'
6364
6503
          
6365
6504
          // Update dmats using an inner product.
6366
6505
          if (combinations[r][s] == 0)
6372
6511
              for (unsigned int tu = 0; tu < 6; tu++)
6373
6512
              {
6374
6513
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
6375
 
              }// end loop over 'tu'
6376
 
            }// end loop over 'u'
6377
 
          }// end loop over 't'
 
6514
              } // end loop over 'tu'
 
6515
            } // end loop over 'u'
 
6516
          } // end loop over 't'
6378
6517
          }
6379
6518
          
6380
6519
          if (combinations[r][s] == 1)
6386
6525
              for (unsigned int tu = 0; tu < 6; tu++)
6387
6526
              {
6388
6527
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
6389
 
              }// end loop over 'tu'
6390
 
            }// end loop over 'u'
6391
 
          }// end loop over 't'
 
6528
              } // end loop over 'tu'
 
6529
            } // end loop over 'u'
 
6530
          } // end loop over 't'
6392
6531
          }
6393
6532
          
6394
 
        }// end loop over 's'
 
6533
        } // end loop over 's'
6395
6534
        for (unsigned int s = 0; s < 6; s++)
6396
6535
        {
6397
6536
          for (unsigned int t = 0; t < 6; t++)
6398
6537
          {
6399
6538
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
6400
 
          }// end loop over 't'
6401
 
        }// end loop over 's'
6402
 
      }// end loop over 'r'
 
6539
          } // end loop over 't'
 
6540
        } // end loop over 's'
 
6541
      } // end loop over 'r'
6403
6542
      
6404
6543
      // Transform derivatives back to physical element
6405
6544
      for (unsigned int r = 0; r < num_derivatives; r++)
6407
6546
        for (unsigned int s = 0; s < num_derivatives; s++)
6408
6547
        {
6409
6548
          values[r] += transform[r][s]*derivatives[s];
6410
 
        }// end loop over 's'
6411
 
      }// end loop over 'r'
 
6549
        } // end loop over 's'
 
6550
      } // end loop over 'r'
6412
6551
        break;
6413
6552
      }
6414
6553
    case 3:
6463
6602
      for (unsigned int r = 0; r < 4; r++)
6464
6603
      {
6465
6604
        derivatives[r] = 0.0;
6466
 
      }// end loop over 'r'
 
6605
      } // end loop over 'r'
6467
6606
      
6468
6607
      // Declare derivative matrix (of polynomial basis).
6469
6608
      double dmats[6][6] = \
6497
6636
            dmats[t][u] = 1.0;
6498
6637
            }
6499
6638
            
6500
 
          }// end loop over 'u'
6501
 
        }// end loop over 't'
 
6639
          } // end loop over 'u'
 
6640
        } // end loop over 't'
6502
6641
        
6503
6642
        // Looping derivative order to generate dmats.
6504
6643
        for (unsigned int s = 0; s < n; s++)
6510
6649
            {
6511
6650
              dmats_old[t][u] = dmats[t][u];
6512
6651
              dmats[t][u] = 0.0;
6513
 
            }// end loop over 'u'
6514
 
          }// end loop over 't'
 
6652
            } // end loop over 'u'
 
6653
          } // end loop over 't'
6515
6654
          
6516
6655
          // Update dmats using an inner product.
6517
6656
          if (combinations[r][s] == 0)
6523
6662
              for (unsigned int tu = 0; tu < 6; tu++)
6524
6663
              {
6525
6664
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
6526
 
              }// end loop over 'tu'
6527
 
            }// end loop over 'u'
6528
 
          }// end loop over 't'
 
6665
              } // end loop over 'tu'
 
6666
            } // end loop over 'u'
 
6667
          } // end loop over 't'
6529
6668
          }
6530
6669
          
6531
6670
          if (combinations[r][s] == 1)
6537
6676
              for (unsigned int tu = 0; tu < 6; tu++)
6538
6677
              {
6539
6678
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
6540
 
              }// end loop over 'tu'
6541
 
            }// end loop over 'u'
6542
 
          }// end loop over 't'
 
6679
              } // end loop over 'tu'
 
6680
            } // end loop over 'u'
 
6681
          } // end loop over 't'
6543
6682
          }
6544
6683
          
6545
 
        }// end loop over 's'
 
6684
        } // end loop over 's'
6546
6685
        for (unsigned int s = 0; s < 6; s++)
6547
6686
        {
6548
6687
          for (unsigned int t = 0; t < 6; t++)
6549
6688
          {
6550
6689
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
6551
 
          }// end loop over 't'
6552
 
        }// end loop over 's'
6553
 
      }// end loop over 'r'
 
6690
          } // end loop over 't'
 
6691
        } // end loop over 's'
 
6692
      } // end loop over 'r'
6554
6693
      
6555
6694
      // Transform derivatives back to physical element
6556
6695
      for (unsigned int r = 0; r < num_derivatives; r++)
6558
6697
        for (unsigned int s = 0; s < num_derivatives; s++)
6559
6698
        {
6560
6699
          values[r] += transform[r][s]*derivatives[s];
6561
 
        }// end loop over 's'
6562
 
      }// end loop over 'r'
 
6700
        } // end loop over 's'
 
6701
      } // end loop over 'r'
6563
6702
        break;
6564
6703
      }
6565
6704
    case 4:
6614
6753
      for (unsigned int r = 0; r < 4; r++)
6615
6754
      {
6616
6755
        derivatives[r] = 0.0;
6617
 
      }// end loop over 'r'
 
6756
      } // end loop over 'r'
6618
6757
      
6619
6758
      // Declare derivative matrix (of polynomial basis).
6620
6759
      double dmats[6][6] = \
6648
6787
            dmats[t][u] = 1.0;
6649
6788
            }
6650
6789
            
6651
 
          }// end loop over 'u'
6652
 
        }// end loop over 't'
 
6790
          } // end loop over 'u'
 
6791
        } // end loop over 't'
6653
6792
        
6654
6793
        // Looping derivative order to generate dmats.
6655
6794
        for (unsigned int s = 0; s < n; s++)
6661
6800
            {
6662
6801
              dmats_old[t][u] = dmats[t][u];
6663
6802
              dmats[t][u] = 0.0;
6664
 
            }// end loop over 'u'
6665
 
          }// end loop over 't'
 
6803
            } // end loop over 'u'
 
6804
          } // end loop over 't'
6666
6805
          
6667
6806
          // Update dmats using an inner product.
6668
6807
          if (combinations[r][s] == 0)
6674
6813
              for (unsigned int tu = 0; tu < 6; tu++)
6675
6814
              {
6676
6815
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
6677
 
              }// end loop over 'tu'
6678
 
            }// end loop over 'u'
6679
 
          }// end loop over 't'
 
6816
              } // end loop over 'tu'
 
6817
            } // end loop over 'u'
 
6818
          } // end loop over 't'
6680
6819
          }
6681
6820
          
6682
6821
          if (combinations[r][s] == 1)
6688
6827
              for (unsigned int tu = 0; tu < 6; tu++)
6689
6828
              {
6690
6829
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
6691
 
              }// end loop over 'tu'
6692
 
            }// end loop over 'u'
6693
 
          }// end loop over 't'
 
6830
              } // end loop over 'tu'
 
6831
            } // end loop over 'u'
 
6832
          } // end loop over 't'
6694
6833
          }
6695
6834
          
6696
 
        }// end loop over 's'
 
6835
        } // end loop over 's'
6697
6836
        for (unsigned int s = 0; s < 6; s++)
6698
6837
        {
6699
6838
          for (unsigned int t = 0; t < 6; t++)
6700
6839
          {
6701
6840
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
6702
 
          }// end loop over 't'
6703
 
        }// end loop over 's'
6704
 
      }// end loop over 'r'
 
6841
          } // end loop over 't'
 
6842
        } // end loop over 's'
 
6843
      } // end loop over 'r'
6705
6844
      
6706
6845
      // Transform derivatives back to physical element
6707
6846
      for (unsigned int r = 0; r < num_derivatives; r++)
6709
6848
        for (unsigned int s = 0; s < num_derivatives; s++)
6710
6849
        {
6711
6850
          values[r] += transform[r][s]*derivatives[s];
6712
 
        }// end loop over 's'
6713
 
      }// end loop over 'r'
 
6851
        } // end loop over 's'
 
6852
      } // end loop over 'r'
6714
6853
        break;
6715
6854
      }
6716
6855
    case 5:
6765
6904
      for (unsigned int r = 0; r < 4; r++)
6766
6905
      {
6767
6906
        derivatives[r] = 0.0;
6768
 
      }// end loop over 'r'
 
6907
      } // end loop over 'r'
6769
6908
      
6770
6909
      // Declare derivative matrix (of polynomial basis).
6771
6910
      double dmats[6][6] = \
6799
6938
            dmats[t][u] = 1.0;
6800
6939
            }
6801
6940
            
6802
 
          }// end loop over 'u'
6803
 
        }// end loop over 't'
 
6941
          } // end loop over 'u'
 
6942
        } // end loop over 't'
6804
6943
        
6805
6944
        // Looping derivative order to generate dmats.
6806
6945
        for (unsigned int s = 0; s < n; s++)
6812
6951
            {
6813
6952
              dmats_old[t][u] = dmats[t][u];
6814
6953
              dmats[t][u] = 0.0;
6815
 
            }// end loop over 'u'
6816
 
          }// end loop over 't'
 
6954
            } // end loop over 'u'
 
6955
          } // end loop over 't'
6817
6956
          
6818
6957
          // Update dmats using an inner product.
6819
6958
          if (combinations[r][s] == 0)
6825
6964
              for (unsigned int tu = 0; tu < 6; tu++)
6826
6965
              {
6827
6966
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
6828
 
              }// end loop over 'tu'
6829
 
            }// end loop over 'u'
6830
 
          }// end loop over 't'
 
6967
              } // end loop over 'tu'
 
6968
            } // end loop over 'u'
 
6969
          } // end loop over 't'
6831
6970
          }
6832
6971
          
6833
6972
          if (combinations[r][s] == 1)
6839
6978
              for (unsigned int tu = 0; tu < 6; tu++)
6840
6979
              {
6841
6980
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
6842
 
              }// end loop over 'tu'
6843
 
            }// end loop over 'u'
6844
 
          }// end loop over 't'
 
6981
              } // end loop over 'tu'
 
6982
            } // end loop over 'u'
 
6983
          } // end loop over 't'
6845
6984
          }
6846
6985
          
6847
 
        }// end loop over 's'
 
6986
        } // end loop over 's'
6848
6987
        for (unsigned int s = 0; s < 6; s++)
6849
6988
        {
6850
6989
          for (unsigned int t = 0; t < 6; t++)
6851
6990
          {
6852
6991
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
6853
 
          }// end loop over 't'
6854
 
        }// end loop over 's'
6855
 
      }// end loop over 'r'
 
6992
          } // end loop over 't'
 
6993
        } // end loop over 's'
 
6994
      } // end loop over 'r'
6856
6995
      
6857
6996
      // Transform derivatives back to physical element
6858
6997
      for (unsigned int r = 0; r < num_derivatives; r++)
6860
6999
        for (unsigned int s = 0; s < num_derivatives; s++)
6861
7000
        {
6862
7001
          values[r] += transform[r][s]*derivatives[s];
6863
 
        }// end loop over 's'
6864
 
      }// end loop over 'r'
 
7002
        } // end loop over 's'
 
7003
      } // end loop over 'r'
6865
7004
        break;
6866
7005
      }
6867
7006
    case 6:
6916
7055
      for (unsigned int r = 0; r < 4; r++)
6917
7056
      {
6918
7057
        derivatives[r] = 0.0;
6919
 
      }// end loop over 'r'
 
7058
      } // end loop over 'r'
6920
7059
      
6921
7060
      // Declare derivative matrix (of polynomial basis).
6922
7061
      double dmats[6][6] = \
6950
7089
            dmats[t][u] = 1.0;
6951
7090
            }
6952
7091
            
6953
 
          }// end loop over 'u'
6954
 
        }// end loop over 't'
 
7092
          } // end loop over 'u'
 
7093
        } // end loop over 't'
6955
7094
        
6956
7095
        // Looping derivative order to generate dmats.
6957
7096
        for (unsigned int s = 0; s < n; s++)
6963
7102
            {
6964
7103
              dmats_old[t][u] = dmats[t][u];
6965
7104
              dmats[t][u] = 0.0;
6966
 
            }// end loop over 'u'
6967
 
          }// end loop over 't'
 
7105
            } // end loop over 'u'
 
7106
          } // end loop over 't'
6968
7107
          
6969
7108
          // Update dmats using an inner product.
6970
7109
          if (combinations[r][s] == 0)
6976
7115
              for (unsigned int tu = 0; tu < 6; tu++)
6977
7116
              {
6978
7117
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
6979
 
              }// end loop over 'tu'
6980
 
            }// end loop over 'u'
6981
 
          }// end loop over 't'
 
7118
              } // end loop over 'tu'
 
7119
            } // end loop over 'u'
 
7120
          } // end loop over 't'
6982
7121
          }
6983
7122
          
6984
7123
          if (combinations[r][s] == 1)
6990
7129
              for (unsigned int tu = 0; tu < 6; tu++)
6991
7130
              {
6992
7131
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
6993
 
              }// end loop over 'tu'
6994
 
            }// end loop over 'u'
6995
 
          }// end loop over 't'
 
7132
              } // end loop over 'tu'
 
7133
            } // end loop over 'u'
 
7134
          } // end loop over 't'
6996
7135
          }
6997
7136
          
6998
 
        }// end loop over 's'
 
7137
        } // end loop over 's'
6999
7138
        for (unsigned int s = 0; s < 6; s++)
7000
7139
        {
7001
7140
          for (unsigned int t = 0; t < 6; t++)
7002
7141
          {
7003
7142
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
7004
 
          }// end loop over 't'
7005
 
        }// end loop over 's'
7006
 
      }// end loop over 'r'
 
7143
          } // end loop over 't'
 
7144
        } // end loop over 's'
 
7145
      } // end loop over 'r'
7007
7146
      
7008
7147
      // Transform derivatives back to physical element
7009
7148
      for (unsigned int r = 0; r < num_derivatives; r++)
7011
7150
        for (unsigned int s = 0; s < num_derivatives; s++)
7012
7151
        {
7013
7152
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
7014
 
        }// end loop over 's'
7015
 
      }// end loop over 'r'
 
7153
        } // end loop over 's'
 
7154
      } // end loop over 'r'
7016
7155
        break;
7017
7156
      }
7018
7157
    case 7:
7067
7206
      for (unsigned int r = 0; r < 4; r++)
7068
7207
      {
7069
7208
        derivatives[r] = 0.0;
7070
 
      }// end loop over 'r'
 
7209
      } // end loop over 'r'
7071
7210
      
7072
7211
      // Declare derivative matrix (of polynomial basis).
7073
7212
      double dmats[6][6] = \
7101
7240
            dmats[t][u] = 1.0;
7102
7241
            }
7103
7242
            
7104
 
          }// end loop over 'u'
7105
 
        }// end loop over 't'
 
7243
          } // end loop over 'u'
 
7244
        } // end loop over 't'
7106
7245
        
7107
7246
        // Looping derivative order to generate dmats.
7108
7247
        for (unsigned int s = 0; s < n; s++)
7114
7253
            {
7115
7254
              dmats_old[t][u] = dmats[t][u];
7116
7255
              dmats[t][u] = 0.0;
7117
 
            }// end loop over 'u'
7118
 
          }// end loop over 't'
 
7256
            } // end loop over 'u'
 
7257
          } // end loop over 't'
7119
7258
          
7120
7259
          // Update dmats using an inner product.
7121
7260
          if (combinations[r][s] == 0)
7127
7266
              for (unsigned int tu = 0; tu < 6; tu++)
7128
7267
              {
7129
7268
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
7130
 
              }// end loop over 'tu'
7131
 
            }// end loop over 'u'
7132
 
          }// end loop over 't'
 
7269
              } // end loop over 'tu'
 
7270
            } // end loop over 'u'
 
7271
          } // end loop over 't'
7133
7272
          }
7134
7273
          
7135
7274
          if (combinations[r][s] == 1)
7141
7280
              for (unsigned int tu = 0; tu < 6; tu++)
7142
7281
              {
7143
7282
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
7144
 
              }// end loop over 'tu'
7145
 
            }// end loop over 'u'
7146
 
          }// end loop over 't'
 
7283
              } // end loop over 'tu'
 
7284
            } // end loop over 'u'
 
7285
          } // end loop over 't'
7147
7286
          }
7148
7287
          
7149
 
        }// end loop over 's'
 
7288
        } // end loop over 's'
7150
7289
        for (unsigned int s = 0; s < 6; s++)
7151
7290
        {
7152
7291
          for (unsigned int t = 0; t < 6; t++)
7153
7292
          {
7154
7293
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
7155
 
          }// end loop over 't'
7156
 
        }// end loop over 's'
7157
 
      }// end loop over 'r'
 
7294
          } // end loop over 't'
 
7295
        } // end loop over 's'
 
7296
      } // end loop over 'r'
7158
7297
      
7159
7298
      // Transform derivatives back to physical element
7160
7299
      for (unsigned int r = 0; r < num_derivatives; r++)
7162
7301
        for (unsigned int s = 0; s < num_derivatives; s++)
7163
7302
        {
7164
7303
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
7165
 
        }// end loop over 's'
7166
 
      }// end loop over 'r'
 
7304
        } // end loop over 's'
 
7305
      } // end loop over 'r'
7167
7306
        break;
7168
7307
      }
7169
7308
    case 8:
7218
7357
      for (unsigned int r = 0; r < 4; r++)
7219
7358
      {
7220
7359
        derivatives[r] = 0.0;
7221
 
      }// end loop over 'r'
 
7360
      } // end loop over 'r'
7222
7361
      
7223
7362
      // Declare derivative matrix (of polynomial basis).
7224
7363
      double dmats[6][6] = \
7252
7391
            dmats[t][u] = 1.0;
7253
7392
            }
7254
7393
            
7255
 
          }// end loop over 'u'
7256
 
        }// end loop over 't'
 
7394
          } // end loop over 'u'
 
7395
        } // end loop over 't'
7257
7396
        
7258
7397
        // Looping derivative order to generate dmats.
7259
7398
        for (unsigned int s = 0; s < n; s++)
7265
7404
            {
7266
7405
              dmats_old[t][u] = dmats[t][u];
7267
7406
              dmats[t][u] = 0.0;
7268
 
            }// end loop over 'u'
7269
 
          }// end loop over 't'
 
7407
            } // end loop over 'u'
 
7408
          } // end loop over 't'
7270
7409
          
7271
7410
          // Update dmats using an inner product.
7272
7411
          if (combinations[r][s] == 0)
7278
7417
              for (unsigned int tu = 0; tu < 6; tu++)
7279
7418
              {
7280
7419
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
7281
 
              }// end loop over 'tu'
7282
 
            }// end loop over 'u'
7283
 
          }// end loop over 't'
 
7420
              } // end loop over 'tu'
 
7421
            } // end loop over 'u'
 
7422
          } // end loop over 't'
7284
7423
          }
7285
7424
          
7286
7425
          if (combinations[r][s] == 1)
7292
7431
              for (unsigned int tu = 0; tu < 6; tu++)
7293
7432
              {
7294
7433
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
7295
 
              }// end loop over 'tu'
7296
 
            }// end loop over 'u'
7297
 
          }// end loop over 't'
 
7434
              } // end loop over 'tu'
 
7435
            } // end loop over 'u'
 
7436
          } // end loop over 't'
7298
7437
          }
7299
7438
          
7300
 
        }// end loop over 's'
 
7439
        } // end loop over 's'
7301
7440
        for (unsigned int s = 0; s < 6; s++)
7302
7441
        {
7303
7442
          for (unsigned int t = 0; t < 6; t++)
7304
7443
          {
7305
7444
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
7306
 
          }// end loop over 't'
7307
 
        }// end loop over 's'
7308
 
      }// end loop over 'r'
 
7445
          } // end loop over 't'
 
7446
        } // end loop over 's'
 
7447
      } // end loop over 'r'
7309
7448
      
7310
7449
      // Transform derivatives back to physical element
7311
7450
      for (unsigned int r = 0; r < num_derivatives; r++)
7313
7452
        for (unsigned int s = 0; s < num_derivatives; s++)
7314
7453
        {
7315
7454
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
7316
 
        }// end loop over 's'
7317
 
      }// end loop over 'r'
 
7455
        } // end loop over 's'
 
7456
      } // end loop over 'r'
7318
7457
        break;
7319
7458
      }
7320
7459
    case 9:
7369
7508
      for (unsigned int r = 0; r < 4; r++)
7370
7509
      {
7371
7510
        derivatives[r] = 0.0;
7372
 
      }// end loop over 'r'
 
7511
      } // end loop over 'r'
7373
7512
      
7374
7513
      // Declare derivative matrix (of polynomial basis).
7375
7514
      double dmats[6][6] = \
7403
7542
            dmats[t][u] = 1.0;
7404
7543
            }
7405
7544
            
7406
 
          }// end loop over 'u'
7407
 
        }// end loop over 't'
 
7545
          } // end loop over 'u'
 
7546
        } // end loop over 't'
7408
7547
        
7409
7548
        // Looping derivative order to generate dmats.
7410
7549
        for (unsigned int s = 0; s < n; s++)
7416
7555
            {
7417
7556
              dmats_old[t][u] = dmats[t][u];
7418
7557
              dmats[t][u] = 0.0;
7419
 
            }// end loop over 'u'
7420
 
          }// end loop over 't'
 
7558
            } // end loop over 'u'
 
7559
          } // end loop over 't'
7421
7560
          
7422
7561
          // Update dmats using an inner product.
7423
7562
          if (combinations[r][s] == 0)
7429
7568
              for (unsigned int tu = 0; tu < 6; tu++)
7430
7569
              {
7431
7570
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
7432
 
              }// end loop over 'tu'
7433
 
            }// end loop over 'u'
7434
 
          }// end loop over 't'
 
7571
              } // end loop over 'tu'
 
7572
            } // end loop over 'u'
 
7573
          } // end loop over 't'
7435
7574
          }
7436
7575
          
7437
7576
          if (combinations[r][s] == 1)
7443
7582
              for (unsigned int tu = 0; tu < 6; tu++)
7444
7583
              {
7445
7584
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
7446
 
              }// end loop over 'tu'
7447
 
            }// end loop over 'u'
7448
 
          }// end loop over 't'
 
7585
              } // end loop over 'tu'
 
7586
            } // end loop over 'u'
 
7587
          } // end loop over 't'
7449
7588
          }
7450
7589
          
7451
 
        }// end loop over 's'
 
7590
        } // end loop over 's'
7452
7591
        for (unsigned int s = 0; s < 6; s++)
7453
7592
        {
7454
7593
          for (unsigned int t = 0; t < 6; t++)
7455
7594
          {
7456
7595
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
7457
 
          }// end loop over 't'
7458
 
        }// end loop over 's'
7459
 
      }// end loop over 'r'
 
7596
          } // end loop over 't'
 
7597
        } // end loop over 's'
 
7598
      } // end loop over 'r'
7460
7599
      
7461
7600
      // Transform derivatives back to physical element
7462
7601
      for (unsigned int r = 0; r < num_derivatives; r++)
7464
7603
        for (unsigned int s = 0; s < num_derivatives; s++)
7465
7604
        {
7466
7605
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
7467
 
        }// end loop over 's'
7468
 
      }// end loop over 'r'
 
7606
        } // end loop over 's'
 
7607
      } // end loop over 'r'
7469
7608
        break;
7470
7609
      }
7471
7610
    case 10:
7520
7659
      for (unsigned int r = 0; r < 4; r++)
7521
7660
      {
7522
7661
        derivatives[r] = 0.0;
7523
 
      }// end loop over 'r'
 
7662
      } // end loop over 'r'
7524
7663
      
7525
7664
      // Declare derivative matrix (of polynomial basis).
7526
7665
      double dmats[6][6] = \
7554
7693
            dmats[t][u] = 1.0;
7555
7694
            }
7556
7695
            
7557
 
          }// end loop over 'u'
7558
 
        }// end loop over 't'
 
7696
          } // end loop over 'u'
 
7697
        } // end loop over 't'
7559
7698
        
7560
7699
        // Looping derivative order to generate dmats.
7561
7700
        for (unsigned int s = 0; s < n; s++)
7567
7706
            {
7568
7707
              dmats_old[t][u] = dmats[t][u];
7569
7708
              dmats[t][u] = 0.0;
7570
 
            }// end loop over 'u'
7571
 
          }// end loop over 't'
 
7709
            } // end loop over 'u'
 
7710
          } // end loop over 't'
7572
7711
          
7573
7712
          // Update dmats using an inner product.
7574
7713
          if (combinations[r][s] == 0)
7580
7719
              for (unsigned int tu = 0; tu < 6; tu++)
7581
7720
              {
7582
7721
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
7583
 
              }// end loop over 'tu'
7584
 
            }// end loop over 'u'
7585
 
          }// end loop over 't'
 
7722
              } // end loop over 'tu'
 
7723
            } // end loop over 'u'
 
7724
          } // end loop over 't'
7586
7725
          }
7587
7726
          
7588
7727
          if (combinations[r][s] == 1)
7594
7733
              for (unsigned int tu = 0; tu < 6; tu++)
7595
7734
              {
7596
7735
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
7597
 
              }// end loop over 'tu'
7598
 
            }// end loop over 'u'
7599
 
          }// end loop over 't'
 
7736
              } // end loop over 'tu'
 
7737
            } // end loop over 'u'
 
7738
          } // end loop over 't'
7600
7739
          }
7601
7740
          
7602
 
        }// end loop over 's'
 
7741
        } // end loop over 's'
7603
7742
        for (unsigned int s = 0; s < 6; s++)
7604
7743
        {
7605
7744
          for (unsigned int t = 0; t < 6; t++)
7606
7745
          {
7607
7746
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
7608
 
          }// end loop over 't'
7609
 
        }// end loop over 's'
7610
 
      }// end loop over 'r'
 
7747
          } // end loop over 't'
 
7748
        } // end loop over 's'
 
7749
      } // end loop over 'r'
7611
7750
      
7612
7751
      // Transform derivatives back to physical element
7613
7752
      for (unsigned int r = 0; r < num_derivatives; r++)
7615
7754
        for (unsigned int s = 0; s < num_derivatives; s++)
7616
7755
        {
7617
7756
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
7618
 
        }// end loop over 's'
7619
 
      }// end loop over 'r'
 
7757
        } // end loop over 's'
 
7758
      } // end loop over 'r'
7620
7759
        break;
7621
7760
      }
7622
7761
    case 11:
7671
7810
      for (unsigned int r = 0; r < 4; r++)
7672
7811
      {
7673
7812
        derivatives[r] = 0.0;
7674
 
      }// end loop over 'r'
 
7813
      } // end loop over 'r'
7675
7814
      
7676
7815
      // Declare derivative matrix (of polynomial basis).
7677
7816
      double dmats[6][6] = \
7705
7844
            dmats[t][u] = 1.0;
7706
7845
            }
7707
7846
            
7708
 
          }// end loop over 'u'
7709
 
        }// end loop over 't'
 
7847
          } // end loop over 'u'
 
7848
        } // end loop over 't'
7710
7849
        
7711
7850
        // Looping derivative order to generate dmats.
7712
7851
        for (unsigned int s = 0; s < n; s++)
7718
7857
            {
7719
7858
              dmats_old[t][u] = dmats[t][u];
7720
7859
              dmats[t][u] = 0.0;
7721
 
            }// end loop over 'u'
7722
 
          }// end loop over 't'
 
7860
            } // end loop over 'u'
 
7861
          } // end loop over 't'
7723
7862
          
7724
7863
          // Update dmats using an inner product.
7725
7864
          if (combinations[r][s] == 0)
7731
7870
              for (unsigned int tu = 0; tu < 6; tu++)
7732
7871
              {
7733
7872
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
7734
 
              }// end loop over 'tu'
7735
 
            }// end loop over 'u'
7736
 
          }// end loop over 't'
 
7873
              } // end loop over 'tu'
 
7874
            } // end loop over 'u'
 
7875
          } // end loop over 't'
7737
7876
          }
7738
7877
          
7739
7878
          if (combinations[r][s] == 1)
7745
7884
              for (unsigned int tu = 0; tu < 6; tu++)
7746
7885
              {
7747
7886
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
7748
 
              }// end loop over 'tu'
7749
 
            }// end loop over 'u'
7750
 
          }// end loop over 't'
 
7887
              } // end loop over 'tu'
 
7888
            } // end loop over 'u'
 
7889
          } // end loop over 't'
7751
7890
          }
7752
7891
          
7753
 
        }// end loop over 's'
 
7892
        } // end loop over 's'
7754
7893
        for (unsigned int s = 0; s < 6; s++)
7755
7894
        {
7756
7895
          for (unsigned int t = 0; t < 6; t++)
7757
7896
          {
7758
7897
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
7759
 
          }// end loop over 't'
7760
 
        }// end loop over 's'
7761
 
      }// end loop over 'r'
 
7898
          } // end loop over 't'
 
7899
        } // end loop over 's'
 
7900
      } // end loop over 'r'
7762
7901
      
7763
7902
      // Transform derivatives back to physical element
7764
7903
      for (unsigned int r = 0; r < num_derivatives; r++)
7766
7905
        for (unsigned int s = 0; s < num_derivatives; s++)
7767
7906
        {
7768
7907
          values[num_derivatives + r] += transform[r][s]*derivatives[s];
7769
 
        }// end loop over 's'
7770
 
      }// end loop over 'r'
 
7908
        } // end loop over 's'
 
7909
      } // end loop over 'r'
7771
7910
        break;
7772
7911
      }
7773
7912
    case 12:
7808
7947
      for (unsigned int r = 0; r < 4; r++)
7809
7948
      {
7810
7949
        derivatives[r] = 0.0;
7811
 
      }// end loop over 'r'
 
7950
      } // end loop over 'r'
7812
7951
      
7813
7952
      // Declare derivative matrix (of polynomial basis).
7814
7953
      double dmats[3][3] = \
7836
7975
            dmats[t][u] = 1.0;
7837
7976
            }
7838
7977
            
7839
 
          }// end loop over 'u'
7840
 
        }// end loop over 't'
 
7978
          } // end loop over 'u'
 
7979
        } // end loop over 't'
7841
7980
        
7842
7981
        // Looping derivative order to generate dmats.
7843
7982
        for (unsigned int s = 0; s < n; s++)
7849
7988
            {
7850
7989
              dmats_old[t][u] = dmats[t][u];
7851
7990
              dmats[t][u] = 0.0;
7852
 
            }// end loop over 'u'
7853
 
          }// end loop over 't'
 
7991
            } // end loop over 'u'
 
7992
          } // end loop over 't'
7854
7993
          
7855
7994
          // Update dmats using an inner product.
7856
7995
          if (combinations[r][s] == 0)
7862
8001
              for (unsigned int tu = 0; tu < 3; tu++)
7863
8002
              {
7864
8003
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
7865
 
              }// end loop over 'tu'
7866
 
            }// end loop over 'u'
7867
 
          }// end loop over 't'
 
8004
              } // end loop over 'tu'
 
8005
            } // end loop over 'u'
 
8006
          } // end loop over 't'
7868
8007
          }
7869
8008
          
7870
8009
          if (combinations[r][s] == 1)
7876
8015
              for (unsigned int tu = 0; tu < 3; tu++)
7877
8016
              {
7878
8017
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
7879
 
              }// end loop over 'tu'
7880
 
            }// end loop over 'u'
7881
 
          }// end loop over 't'
 
8018
              } // end loop over 'tu'
 
8019
            } // end loop over 'u'
 
8020
          } // end loop over 't'
7882
8021
          }
7883
8022
          
7884
 
        }// end loop over 's'
 
8023
        } // end loop over 's'
7885
8024
        for (unsigned int s = 0; s < 3; s++)
7886
8025
        {
7887
8026
          for (unsigned int t = 0; t < 3; t++)
7888
8027
          {
7889
8028
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
7890
 
          }// end loop over 't'
7891
 
        }// end loop over 's'
7892
 
      }// end loop over 'r'
 
8029
          } // end loop over 't'
 
8030
        } // end loop over 's'
 
8031
      } // end loop over 'r'
7893
8032
      
7894
8033
      // Transform derivatives back to physical element
7895
8034
      for (unsigned int r = 0; r < num_derivatives; r++)
7897
8036
        for (unsigned int s = 0; s < num_derivatives; s++)
7898
8037
        {
7899
8038
          values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
7900
 
        }// end loop over 's'
7901
 
      }// end loop over 'r'
 
8039
        } // end loop over 's'
 
8040
      } // end loop over 'r'
7902
8041
        break;
7903
8042
      }
7904
8043
    case 13:
7939
8078
      for (unsigned int r = 0; r < 4; r++)
7940
8079
      {
7941
8080
        derivatives[r] = 0.0;
7942
 
      }// end loop over 'r'
 
8081
      } // end loop over 'r'
7943
8082
      
7944
8083
      // Declare derivative matrix (of polynomial basis).
7945
8084
      double dmats[3][3] = \
7967
8106
            dmats[t][u] = 1.0;
7968
8107
            }
7969
8108
            
7970
 
          }// end loop over 'u'
7971
 
        }// end loop over 't'
 
8109
          } // end loop over 'u'
 
8110
        } // end loop over 't'
7972
8111
        
7973
8112
        // Looping derivative order to generate dmats.
7974
8113
        for (unsigned int s = 0; s < n; s++)
7980
8119
            {
7981
8120
              dmats_old[t][u] = dmats[t][u];
7982
8121
              dmats[t][u] = 0.0;
7983
 
            }// end loop over 'u'
7984
 
          }// end loop over 't'
 
8122
            } // end loop over 'u'
 
8123
          } // end loop over 't'
7985
8124
          
7986
8125
          // Update dmats using an inner product.
7987
8126
          if (combinations[r][s] == 0)
7993
8132
              for (unsigned int tu = 0; tu < 3; tu++)
7994
8133
              {
7995
8134
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
7996
 
              }// end loop over 'tu'
7997
 
            }// end loop over 'u'
7998
 
          }// end loop over 't'
 
8135
              } // end loop over 'tu'
 
8136
            } // end loop over 'u'
 
8137
          } // end loop over 't'
7999
8138
          }
8000
8139
          
8001
8140
          if (combinations[r][s] == 1)
8007
8146
              for (unsigned int tu = 0; tu < 3; tu++)
8008
8147
              {
8009
8148
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
8010
 
              }// end loop over 'tu'
8011
 
            }// end loop over 'u'
8012
 
          }// end loop over 't'
 
8149
              } // end loop over 'tu'
 
8150
            } // end loop over 'u'
 
8151
          } // end loop over 't'
8013
8152
          }
8014
8153
          
8015
 
        }// end loop over 's'
 
8154
        } // end loop over 's'
8016
8155
        for (unsigned int s = 0; s < 3; s++)
8017
8156
        {
8018
8157
          for (unsigned int t = 0; t < 3; t++)
8019
8158
          {
8020
8159
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
8021
 
          }// end loop over 't'
8022
 
        }// end loop over 's'
8023
 
      }// end loop over 'r'
 
8160
          } // end loop over 't'
 
8161
        } // end loop over 's'
 
8162
      } // end loop over 'r'
8024
8163
      
8025
8164
      // Transform derivatives back to physical element
8026
8165
      for (unsigned int r = 0; r < num_derivatives; r++)
8028
8167
        for (unsigned int s = 0; s < num_derivatives; s++)
8029
8168
        {
8030
8169
          values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
8031
 
        }// end loop over 's'
8032
 
      }// end loop over 'r'
 
8170
        } // end loop over 's'
 
8171
      } // end loop over 'r'
8033
8172
        break;
8034
8173
      }
8035
8174
    case 14:
8070
8209
      for (unsigned int r = 0; r < 4; r++)
8071
8210
      {
8072
8211
        derivatives[r] = 0.0;
8073
 
      }// end loop over 'r'
 
8212
      } // end loop over 'r'
8074
8213
      
8075
8214
      // Declare derivative matrix (of polynomial basis).
8076
8215
      double dmats[3][3] = \
8098
8237
            dmats[t][u] = 1.0;
8099
8238
            }
8100
8239
            
8101
 
          }// end loop over 'u'
8102
 
        }// end loop over 't'
 
8240
          } // end loop over 'u'
 
8241
        } // end loop over 't'
8103
8242
        
8104
8243
        // Looping derivative order to generate dmats.
8105
8244
        for (unsigned int s = 0; s < n; s++)
8111
8250
            {
8112
8251
              dmats_old[t][u] = dmats[t][u];
8113
8252
              dmats[t][u] = 0.0;
8114
 
            }// end loop over 'u'
8115
 
          }// end loop over 't'
 
8253
            } // end loop over 'u'
 
8254
          } // end loop over 't'
8116
8255
          
8117
8256
          // Update dmats using an inner product.
8118
8257
          if (combinations[r][s] == 0)
8124
8263
              for (unsigned int tu = 0; tu < 3; tu++)
8125
8264
              {
8126
8265
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
8127
 
              }// end loop over 'tu'
8128
 
            }// end loop over 'u'
8129
 
          }// end loop over 't'
 
8266
              } // end loop over 'tu'
 
8267
            } // end loop over 'u'
 
8268
          } // end loop over 't'
8130
8269
          }
8131
8270
          
8132
8271
          if (combinations[r][s] == 1)
8138
8277
              for (unsigned int tu = 0; tu < 3; tu++)
8139
8278
              {
8140
8279
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
8141
 
              }// end loop over 'tu'
8142
 
            }// end loop over 'u'
8143
 
          }// end loop over 't'
 
8280
              } // end loop over 'tu'
 
8281
            } // end loop over 'u'
 
8282
          } // end loop over 't'
8144
8283
          }
8145
8284
          
8146
 
        }// end loop over 's'
 
8285
        } // end loop over 's'
8147
8286
        for (unsigned int s = 0; s < 3; s++)
8148
8287
        {
8149
8288
          for (unsigned int t = 0; t < 3; t++)
8150
8289
          {
8151
8290
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
8152
 
          }// end loop over 't'
8153
 
        }// end loop over 's'
8154
 
      }// end loop over 'r'
 
8291
          } // end loop over 't'
 
8292
        } // end loop over 's'
 
8293
      } // end loop over 'r'
8155
8294
      
8156
8295
      // Transform derivatives back to physical element
8157
8296
      for (unsigned int r = 0; r < num_derivatives; r++)
8159
8298
        for (unsigned int s = 0; s < num_derivatives; s++)
8160
8299
        {
8161
8300
          values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
8162
 
        }// end loop over 's'
8163
 
      }// end loop over 'r'
 
8301
        } // end loop over 's'
 
8302
      } // end loop over 'r'
8164
8303
        break;
8165
8304
      }
8166
8305
    }
8167
8306
    
8168
8307
  }
8169
8308
 
8170
 
  /// Evaluate order n derivatives of all basis functions at given point x in cell
8171
 
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
8309
  /// Evaluate order n derivatives of basis function i at given point x in cell (non-static member function)
 
8310
  virtual void evaluate_basis_derivatives(std::size_t i,
 
8311
                                          std::size_t n,
 
8312
                                          double* values,
 
8313
                                          const double* x,
 
8314
                                          const double* vertex_coordinates,
 
8315
                                          int cell_orientation) const
 
8316
  {
 
8317
    _evaluate_basis_derivatives(i, n, values, x, vertex_coordinates, cell_orientation);
 
8318
  }
 
8319
 
 
8320
  /// Evaluate order n derivatives of all basis functions at given point x in cell (actual implementation)
 
8321
  static void _evaluate_basis_derivatives_all(std::size_t n,
8172
8322
                                              double* values,
8173
8323
                                              const double* x,
8174
8324
                                              const double* vertex_coordinates,
8175
 
                                              int cell_orientation) const
 
8325
                                              int cell_orientation)
8176
8326
  {
8177
8327
    // Call evaluate_basis_all if order of derivatives is equal to zero.
8178
8328
    if (n == 0)
8179
8329
    {
8180
 
      evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
 
8330
      _evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
8181
8331
      return ;
8182
8332
    }
8183
8333
    
8186
8336
    for (unsigned int r = 0; r < n; r++)
8187
8337
    {
8188
8338
      num_derivatives *= 2;
8189
 
    }// end loop over 'r'
 
8339
    } // end loop over 'r'
8190
8340
    
8191
8341
    // Set values equal to zero.
8192
8342
    for (unsigned int r = 0; r < 15; r++)
8194
8344
      for (unsigned int s = 0; s < 3*num_derivatives; s++)
8195
8345
      {
8196
8346
        values[r*3*num_derivatives + s] = 0.0;
8197
 
      }// end loop over 's'
8198
 
    }// end loop over 'r'
 
8347
      } // end loop over 's'
 
8348
    } // end loop over 'r'
8199
8349
    
8200
8350
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
8201
8351
    if (n > 2)
8208
8358
    for (unsigned int r = 0; r < 12; r++)
8209
8359
    {
8210
8360
      dof_values[r] = 0.0;
8211
 
    }// end loop over 'r'
 
8361
    } // end loop over 'r'
8212
8362
    
8213
8363
    // Loop dofs and call evaluate_basis_derivatives.
8214
8364
    for (unsigned int r = 0; r < 15; r++)
8215
8365
    {
8216
 
      evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
 
8366
      _evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
8217
8367
      for (unsigned int s = 0; s < 3*num_derivatives; s++)
8218
8368
      {
8219
8369
        values[r*3*num_derivatives + s] = dof_values[s];
8220
 
      }// end loop over 's'
8221
 
    }// end loop over 'r'
 
8370
      } // end loop over 's'
 
8371
    } // end loop over 'r'
 
8372
  }
 
8373
 
 
8374
  /// Evaluate order n derivatives of all basis functions at given point x in cell (non-static member function)
 
8375
  virtual void evaluate_basis_derivatives_all(std::size_t n,
 
8376
                                              double* values,
 
8377
                                              const double* x,
 
8378
                                              const double* vertex_coordinates,
 
8379
                                              int cell_orientation) const
 
8380
  {
 
8381
    _evaluate_basis_derivatives_all(n, values, x, vertex_coordinates, cell_orientation);
8222
8382
  }
8223
8383
 
8224
8384
  /// Evaluate linear functional for dof i on the function f
8527
8687
  /// Return a string identifying the dofmap
8528
8688
  virtual const char* signature() const
8529
8689
  {
8530
 
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, None)";
 
8690
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 2, None)";
8531
8691
  }
8532
8692
 
8533
8693
  /// Return true iff mesh entities of topological dimension d are needed
8794
8954
  /// Return a string identifying the dofmap
8795
8955
  virtual const char* signature() const
8796
8956
  {
8797
 
    return "FFC dofmap for VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None)";
 
8957
    return "FFC dofmap for VectorElement('Lagrange', Domain(Cell('triangle', 2)), 2, 2, None)";
8798
8958
  }
8799
8959
 
8800
8960
  /// Return true iff mesh entities of topological dimension d are needed
9110
9270
  /// Return a string identifying the dofmap
9111
9271
  virtual const char* signature() const
9112
9272
  {
9113
 
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
 
9273
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 1, None)";
9114
9274
  }
9115
9275
 
9116
9276
  /// Return true iff mesh entities of topological dimension d are needed
9339
9499
  /// Return a string identifying the dofmap
9340
9500
  virtual const char* signature() const
9341
9501
  {
9342
 
    return "FFC dofmap for MixedElement(*[VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 2, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (3,) })";
 
9502
    return "FFC dofmap for MixedElement(VectorElement('Lagrange', Domain(Cell('triangle', 2)), 2, 2, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 1, None), **{'value_shape': (3,) })";
9343
9503
  }
9344
9504
 
9345
9505
  /// Return true iff mesh entities of topological dimension d are needed
9672
9832
    // Do nothing
9673
9833
  }
9674
9834
 
 
9835
  /// Tabulate which form coefficients are used by this integral
 
9836
  virtual const std::vector<bool> & enabled_coefficients() const
 
9837
  {
 
9838
    static const std::vector<bool> enabled({});
 
9839
    return enabled;
 
9840
  }
 
9841
 
9675
9842
  /// Tabulate the tensor for the contribution from a local cell
9676
9843
  virtual void tabulate_tensor(double*  A,
9677
9844
                               const double * const *  w,
9985
10152
  /// Return a string identifying the form
9986
10153
  virtual const char* signature() const
9987
10154
  {
9988
 
    return "6805f96321c32af0d26b4a8da27c4997a99275496416f0c5a4bc4b8c0fdda5290a36e29d25b765590a42ab14a8d48072f5c30ede6de761e9292bdf37ffb0f17e";
 
10155
    return "c9e96217bdf0b022e1a9d4f90d63fa9416a1eee23b1a54c975139f2f6b9c7cdd37541ea9b58f1cb326bf42f666b50c6dd5e092cd12a6ee87d8ae82cf4a165465";
9989
10156
  }
9990
10157
 
9991
10158
  /// Return the rank of the global tensor (r)
10024
10191
    return 0;
10025
10192
  }
10026
10193
 
 
10194
  /// Return the number of custom domains
 
10195
  virtual std::size_t num_custom_domains() const
 
10196
  {
 
10197
    return 0;
 
10198
  }
 
10199
 
10027
10200
  /// Return whether the form has any cell integrals
10028
10201
  virtual bool has_cell_integrals() const
10029
10202
  {
10048
10221
    return false;
10049
10222
  }
10050
10223
 
 
10224
  /// Return whether the form has any custom integrals
 
10225
  virtual bool has_custom_integrals() const
 
10226
  {
 
10227
    return false;
 
10228
  }
 
10229
 
10051
10230
  /// Create a new finite element for argument function i
10052
10231
  virtual ufc::finite_element* create_finite_element(std::size_t i) const
10053
10232
  {
10112
10291
    return 0;
10113
10292
  }
10114
10293
 
 
10294
  /// Create a new custom integral on sub domain i
 
10295
  virtual ufc::custom_integral* create_custom_integral(std::size_t i) const
 
10296
  {
 
10297
    return 0;
 
10298
  }
 
10299
 
10115
10300
  /// Create a new cell integral on everywhere else
10116
10301
  virtual ufc::cell_integral* create_default_cell_integral() const
10117
10302
  {
10136
10321
    return 0;
10137
10322
  }
10138
10323
 
 
10324
  /// Create a new custom integral on everywhere else
 
10325
  virtual ufc::custom_integral* create_default_custom_integral() const
 
10326
  {
 
10327
    return 0;
 
10328
  }
 
10329
 
10139
10330
};
10140
10331
 
10141
10332
// DOLFIN wrappers
10167
10358
  // Create standard function space (reference version)
10168
10359
  Form_a_FunctionSpace_0(const dolfin::Mesh& mesh):
10169
10360
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
10170
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10171
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh)))
 
10361
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10362
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh)))
10172
10363
  {
10173
10364
    // Do nothing
10174
10365
  }
10175
10366
 
10176
10367
  // Create standard function space (shared pointer version)
10177
 
  Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
 
10368
  Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh):
10178
10369
    dolfin::FunctionSpace(mesh,
10179
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10180
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh)))
 
10370
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10371
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh)))
10181
10372
  {
10182
10373
    // Do nothing
10183
10374
  }
10187
10378
  // Create standard function space (reference version)
10188
10379
  Form_a_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
10189
10380
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
10190
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10191
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh,
 
10381
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10382
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh,
10192
10383
                              dolfin::reference_to_no_delete_pointer(constrained_domain))))
10193
10384
  {
10194
10385
    // Do nothing
10195
10386
  }
10196
10387
 
10197
10388
  // Create standard function space (shared pointer version)
10198
 
  Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
 
10389
  Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
10199
10390
    dolfin::FunctionSpace(mesh,
10200
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10201
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh, constrained_domain)))
 
10391
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10392
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh, constrained_domain)))
10202
10393
  {
10203
10394
    // Do nothing
10204
10395
  }
10208
10399
  // Create restricted function space (reference version)
10209
10400
  Form_a_FunctionSpace_0(const dolfin::Restriction& restriction):
10210
10401
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
10211
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10212
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
 
10402
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10403
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
10213
10404
                                                                                     reference_to_no_delete_pointer(restriction))))
10214
10405
  {
10215
10406
    // Do nothing
10216
10407
  }
10217
10408
 
10218
10409
  // Create restricted function space (shared pointer version)
10219
 
  Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Restriction> restriction):
 
10410
  Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Restriction> restriction):
10220
10411
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
10221
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10222
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
 
10412
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10413
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
10223
10414
                                                                                     restriction)))
10224
10415
  {
10225
10416
    // Do nothing
10241
10432
  // Create standard function space (reference version)
10242
10433
  Form_a_FunctionSpace_1(const dolfin::Mesh& mesh):
10243
10434
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
10244
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10245
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh)))
 
10435
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10436
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh)))
10246
10437
  {
10247
10438
    // Do nothing
10248
10439
  }
10249
10440
 
10250
10441
  // Create standard function space (shared pointer version)
10251
 
  Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh):
 
10442
  Form_a_FunctionSpace_1(std::shared_ptr<const dolfin::Mesh> mesh):
10252
10443
    dolfin::FunctionSpace(mesh,
10253
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10254
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh)))
 
10444
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10445
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh)))
10255
10446
  {
10256
10447
    // Do nothing
10257
10448
  }
10261
10452
  // Create standard function space (reference version)
10262
10453
  Form_a_FunctionSpace_1(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
10263
10454
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
10264
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10265
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh,
 
10455
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10456
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), mesh,
10266
10457
                              dolfin::reference_to_no_delete_pointer(constrained_domain))))
10267
10458
  {
10268
10459
    // Do nothing
10269
10460
  }
10270
10461
 
10271
10462
  // Create standard function space (shared pointer version)
10272
 
  Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
 
10463
  Form_a_FunctionSpace_1(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
10273
10464
    dolfin::FunctionSpace(mesh,
10274
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10275
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh, constrained_domain)))
 
10465
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10466
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()), *mesh, constrained_domain)))
10276
10467
  {
10277
10468
    // Do nothing
10278
10469
  }
10282
10473
  // Create restricted function space (reference version)
10283
10474
  Form_a_FunctionSpace_1(const dolfin::Restriction& restriction):
10284
10475
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
10285
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10286
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
 
10476
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10477
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
10287
10478
                                                                                     reference_to_no_delete_pointer(restriction))))
10288
10479
  {
10289
10480
    // Do nothing
10290
10481
  }
10291
10482
 
10292
10483
  // Create restricted function space (shared pointer version)
10293
 
  Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Restriction> restriction):
 
10484
  Form_a_FunctionSpace_1(std::shared_ptr<const dolfin::Restriction> restriction):
10294
10485
    dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
10295
 
                          boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
10296
 
                          boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
 
10486
                          std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new thstokes2d_finite_element_3()))),
 
10487
                          std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new thstokes2d_dofmap_3()),
10297
10488
                                                                                     restriction)))
10298
10489
  {
10299
10490
    // Do nothing
10317
10508
    _function_spaces[0] = reference_to_no_delete_pointer(V0);
10318
10509
    _function_spaces[1] = reference_to_no_delete_pointer(V1);
10319
10510
 
10320
 
    _ufc_form = boost::shared_ptr<const ufc::form>(new thstokes2d_form_0());
 
10511
    _ufc_form = std::shared_ptr<const ufc::form>(new thstokes2d_form_0());
10321
10512
  }
10322
10513
 
10323
10514
  // Constructor
10324
 
  Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0):
 
10515
  Form_a(std::shared_ptr<const dolfin::FunctionSpace> V1, std::shared_ptr<const dolfin::FunctionSpace> V0):
10325
10516
    dolfin::Form(2, 0)
10326
10517
  {
10327
10518
    _function_spaces[0] = V0;
10328
10519
    _function_spaces[1] = V1;
10329
10520
 
10330
 
    _ufc_form = boost::shared_ptr<const ufc::form>(new thstokes2d_form_0());
 
10521
    _ufc_form = std::shared_ptr<const ufc::form>(new thstokes2d_form_0());
10331
10522
  }
10332
10523
 
10333
10524
  // Destructor