← Back to branch summary

~ubuntu-branches/ubuntu/karmic/gmsh/karmic

« back to all changes in this revision

Viewing changes to Numeric/FunctionSpace.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Daniel Leidert
  • Date: 2008-05-18 12:46:05 UTC
  • mfrom: (1.2.6 upstream)
  • Revision ID: james.westby@ubuntu.com-20080518124605-716xqbqeo07o497k
Tags: 2.2.0-2
[Christophe Prud'homme]
* Bug fix: "gmsh ships no .desktop", thanks to Vassilis Pandis (Closes:
  #375770). Applied the Ubuntu patch.

[Daniel Leidert]
* debian/control (Vcs-Svn): Fixed.
  (Build-Depends): Use texlive instead of tetex-bin.
* debian/gmsh.doc-base (Section): Fixed accordingly to doc-base (>= 0.8.10).
* debian/rules: Removed some variable declarations, that lead to double
  configuration and seem to be useless.
  (build/gmsh): Try to avoid multiple runs by using a stamp.
  (orig-tarball): Renamed to get-orig-source and changed to use uscan.
* debian/watch: Added.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Numeric/FunctionSpace.cpp
 
1
// $Id: FunctionSpace.cpp,v 1.4 2008-03-20 11:44:09 geuzaine Exp $
 
2
//
 
3
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 
4
//
 
5
// This program is free software; you can redistribute it and/or modify
 
6
// it under the terms of the GNU General Public License as published by
 
7
// the Free Software Foundation; either version 2 of the License, or
 
8
// (at your option) any later version.
 
9
//
 
10
// This program is distributed in the hope that it will be useful,
 
11
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
// GNU General Public License for more details.
 
14
//
 
15
// You should have received a copy of the GNU General Public License
 
16
// along with this program; if not, write to the Free Software
 
17
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 
18
// USA.
 
19
// 
 
20
// Please report all bugs and problems to <gmsh@geuz.org>.
 
21
 
 
22
#include <cmath>
 
23
#include "FunctionSpace.h"
 
24
#include "GmshDefines.h"
 
25
 
 
26
Double_Matrix generatePascalTriangle(int order)
 
27
{
 
28
  Double_Matrix monomials((order + 1) * (order + 2) / 2, 2);
 
29
  int index = 0;
 
30
  for (int i = 0; i <= order; i++) {
 
31
    for (int j = 0; j <= i; j++) {
 
32
      monomials(index, 0) = i - j;
 
33
      monomials(index, 1) = j;
 
34
      index++;
 
35
    }
 
36
  }
 
37
  return monomials;
 
38
}
 
39
 
 
40
// generate the exterior hull of the Pascal triangle for Serendipity element
 
41
 
 
42
Double_Matrix generatePascalSerendipityTriangle(int order)
 
43
{
 
44
  Double_Matrix monomials(3 * order, 2);
 
45
 
 
46
  monomials(0, 0) = 0;
 
47
  monomials(0, 1) = 0;
 
48
  
 
49
  int index = 1;
 
50
  for (int i = 1; i <= order; i++) {
 
51
    if (i == order) {
 
52
      for (int j = 0; j <= i; j++) {
 
53
        monomials(index, 0) = i - j;
 
54
        monomials(index, 1) = j;
 
55
        index++;
 
56
      }
 
57
    }
 
58
    else {
 
59
      monomials(index, 0) = i;
 
60
      monomials(index, 1) = 0;
 
61
      index++;
 
62
      monomials(index, 0) = 0;
 
63
      monomials(index, 1) = i;
 
64
      index++;
 
65
    }
 
66
  }
 
67
  return monomials;
 
68
}
 
69
 
 
70
// generate the monomials subspace of all monomials of order exactly == p
 
71
 
 
72
Double_Matrix generateMonomialSubspace(int dim, int p)
 
73
{
 
74
  Double_Matrix monomials;
 
75
  
 
76
  switch (dim) {
 
77
  case 1:
 
78
    monomials = Double_Matrix(1, 1);
 
79
    monomials(0, 0) = p;
 
80
    break;
 
81
  case 2:
 
82
    monomials = Double_Matrix(p + 1, 2);
 
83
    for (int k = 0; k <= p; k++) {
 
84
      monomials(k, 0) = p - k;
 
85
      monomials(k, 1) = k;
 
86
    }
 
87
    break;
 
88
  case 3:
 
89
    monomials = Double_Matrix((p + 1) * (p + 2) / 2, 3);
 
90
    int index = 0;
 
91
    for (int i = 0; i <= p; i++) {
 
92
      for (int k = 0; k <= p - i; k++) {
 
93
        monomials(index, 0) = p - i - k;
 
94
        monomials(index, 1) = k;
 
95
        monomials(index, 2) = i;
 
96
        index++;
 
97
      }
 
98
    }
 
99
    break;
 
100
  }
 
101
  return monomials;
 
102
}
 
103
 
 
104
// generate external hull of the Pascal tetrahedron
 
105
 
 
106
Double_Matrix generatePascalSerendipityTetrahedron(int order)
 
107
{
 
108
  int nbMonomials = 4 + 6 * std::max(0, order - 1) + 
 
109
    4 * std::max(0, (order - 2) * (order - 1) / 2);
 
110
  Double_Matrix monomials(nbMonomials, 3);
 
111
 
 
112
  // order 0
 
113
  
 
114
  monomials.set_all(0);
 
115
 
 
116
  monomials(0, 0) = 0;
 
117
  monomials(0, 1) = 0;
 
118
  monomials(0, 2) = 0;
 
119
 
 
120
  int index = 1;
 
121
  for (int p = 1; p < order; p++) {
 
122
    for (int i = 0; i < 3; i++) {
 
123
      int j = (i + 1) % 3;
 
124
      int k = (i + 2) % 3;
 
125
      for (int ii = 0; ii < p; ii++) {
 
126
        monomials(index, i) = p - ii;
 
127
        monomials(index, j) = ii;
 
128
        monomials(index, k) = 0;
 
129
        index++;
 
130
      }
 
131
    }
 
132
  }
 
133
  Double_Matrix monomialsMaxOrder = generateMonomialSubspace(3, order);
 
134
  int nbMaxOrder = monomialsMaxOrder.size1();
 
135
    
 
136
  Double_Matrix(monomials.touchSubmatrix(index, nbMaxOrder, 0, 3)).memcpy(monomialsMaxOrder);
 
137
  return monomials;
 
138
}
 
139
 
 
140
// generate Pascal tetrahedron
 
141
 
 
142
Double_Matrix generatePascalTetrahedron(int order)
 
143
{
 
144
  int nbMonomials = (order + 1) * (order + 2) * (order + 3) / 6;
 
145
 
 
146
  Double_Matrix monomials(nbMonomials, 3);
 
147
 
 
148
  int index = 0;
 
149
  for (int p = 0; p <= order; p++) {
 
150
    Double_Matrix monOrder = generateMonomialSubspace(3, p);
 
151
    int nb = monOrder.size1();
 
152
    Double_Matrix(monomials.touchSubmatrix(index, nb, 0, 3)).memcpy(monOrder);
 
153
    index += nb;
 
154
  }
 
155
 
 
156
  return monomials;
 
157
}  
 
158
 
 
159
int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
 
160
int nbdoftriangleserendip(int order) { return 3 * order; }
 
161
 
 
162
void nodepositionface0(int order, double *u,  double *v,  double *w)
 
163
{
 
164
  int ndofT = nbdoftriangle(order);
 
165
  if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
 
166
  
 
167
  u[0]= 0.;    v[0]= 0.;    w[0] = 0.;
 
168
  u[1]= order; v[1]= 0;     w[1] = 0.;
 
169
  u[2]= 0.;    v[2]= order; w[2] = 0.;
 
170
 
 
171
  // edges
 
172
  for (int k = 0; k < (order - 1); k++){
 
173
    u[3 + k] = k + 1; 
 
174
    v[3 + k] =0.; 
 
175
    w[3 + k] = 0.;
 
176
 
 
177
    u[3 + order - 1 + k] = order - 1 - k ; 
 
178
    v[3 + order - 1 + k] = k + 1;
 
179
    w[3 + order - 1 + k] = 0.;
 
180
 
 
181
    u[3 + 2 * (order - 1) + k] = 0.;
 
182
    v[3 + 2 * (order - 1) + k] = order - 1 - k;
 
183
    w[3 + 2 * (order - 1) + k] = 0.;
 
184
  }
 
185
  if (order > 2){
 
186
    int nbdoftemp = nbdoftriangle(order - 3);
 
187
    nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)], 
 
188
                      &w[3 + 3* (order - 1)]);
 
189
    for (int k = 0; k < nbdoftemp; k++){
 
190
      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
191
      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
192
      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
 
193
    }
 
194
  }
 
195
  for (int k = 0; k < ndofT; k++){
 
196
    u[k] = u[k] / order;
 
197
    v[k] = v[k] / order;
 
198
    w[k] = w[k] / order;
 
199
  }
 
200
}
 
201
 
 
202
void nodepositionface1(int order,  double *u,  double *v,  double *w)
 
203
{
 
204
   int ndofT = nbdoftriangle(order);
 
205
   if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
 
206
   
 
207
   u[0]= 0.;    v[0]= 0.;  w[0] = 0.;
 
208
   u[1]= order; v[1]= 0;   w[1] = 0.;
 
209
   u[2]= 0.;    v[2]= 0.;  w[2] = order;
 
210
   // edges
 
211
   for (int k = 0; k < (order - 1); k++){
 
212
     u[3 + k] = k + 1; 
 
213
     v[3 + k] = 0.; 
 
214
     w[3 + k] = 0.;
 
215
     
 
216
     u[3 + order - 1 + k] = order - 1 - k;
 
217
     v[3 + order - 1 + k] = 0.;
 
218
     w[3 + order - 1+ k ] = k + 1; 
 
219
     
 
220
     u[3 + 2 * (order - 1) + k] = 0. ;
 
221
     v[3 + 2 * (order - 1) + k] = 0.;
 
222
     w[3 + 2 * (order - 1) + k] = order - 1 - k; 
 
223
   }
 
224
   if (order > 2){
 
225
     int nbdoftemp = nbdoftriangle(order - 3);
 
226
     nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
 
227
                       &w[3 + 3 * (order - 1)]);
 
228
     for (int k = 0; k < nbdoftemp; k++){
 
229
       u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
230
       v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
 
231
       w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
232
     }
 
233
   }
 
234
   for (int k = 0; k < ndofT; k++){
 
235
     u[k] = u[k] / order;
 
236
     v[k] = v[k] / order;
 
237
     w[k] = w[k] / order;
 
238
   }
 
239
}
 
240
 
 
241
void nodepositionface2(int order,  double *u,  double *v,  double *w)
 
242
{
 
243
   int ndofT = nbdoftriangle(order);
 
244
   if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
 
245
   
 
246
   u[0]= order; v[0]= 0.;    w[0] = 0.;
 
247
   u[1]= 0.;    v[1]= order; w[1] = 0.;
 
248
   u[2]= 0.;    v[2]= 0.;    w[2] = order;
 
249
   // edges
 
250
   for (int k = 0; k < (order - 1); k++){
 
251
     u[3 + k] = order - 1 - k;
 
252
     v[3 + k] = k + 1;
 
253
     w[3 + k] = 0.;
 
254
 
 
255
     u[3 + order - 1 + k] = 0.;
 
256
     v[3 + order - 1 + k] = order - 1 - k;
 
257
     w[3 + order - 1 + k] = k + 1; 
 
258
     
 
259
     u[3 + 2 * (order - 1) + k] = k + 1;
 
260
     v[3 + 2 * (order - 1) + k] = 0.;
 
261
     w[3 + 2 * (order - 1) + k] = order - 1 - k; 
 
262
   }
 
263
   if (order > 2){
 
264
     int nbdoftemp = nbdoftriangle(order - 3);
 
265
     nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
 
266
                       &w[3 + 3 * (order - 1)]);
 
267
     for (int k = 0; k < nbdoftemp; k++){
 
268
       u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
269
       v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
270
       w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
271
     }
 
272
   }
 
273
   for (int k = 0; k < ndofT; k++){
 
274
     u[k] = u[k] / order;
 
275
     v[k] = v[k] / order;
 
276
     w[k] = w[k] / order;
 
277
   }
 
278
}
 
279
 
 
280
void nodepositionface3(int order, double *u, double *v,  double *w)
 
281
{
 
282
   int ndofT = nbdoftriangle(order);
 
283
   if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
 
284
   
 
285
   u[0]= 0.; v[0]= 0.;    w[0] = 0.;
 
286
   u[1]= 0.; v[1]= order; w[1] = 0.;
 
287
   u[2]= 0.; v[2]= 0.;    w[2] = order;
 
288
   // edges
 
289
   for (int k = 0; k < (order - 1); k++){
 
290
     u[3 + k] = 0.;
 
291
     v[3 + k]= k + 1;
 
292
     w[3 + k] = 0.;
 
293
 
 
294
     u[3 + order - 1 + k] = 0.;
 
295
     v[3 + order - 1 + k] = order - 1 - k;
 
296
     w[3 + order - 1 + k] = k + 1; 
 
297
     
 
298
     u[3 + 2 * (order - 1) + k] = 0.;
 
299
     v[3 + 2 * (order - 1) + k] = 0.;
 
300
     w[3 + 2 * (order - 1) + k] = order - 1 - k;
 
301
   }
 
302
   if (order > 2){
 
303
     int nbdoftemp = nbdoftriangle(order - 3);
 
304
     nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
 
305
                       &w[3 + 3 * (order - 1)]);
 
306
     for (int k = 0; k < nbdoftemp; k++){
 
307
       u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
 
308
       v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
309
       w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
 
310
     }
 
311
   }
 
312
   for (int k = 0; k < ndofT; k++){
 
313
     u[k] = u[k] / order;
 
314
     v[k] = v[k] / order;
 
315
     w[k] = w[k] / order;
 
316
   }
 
317
}
 
318
 
 
319
Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip) 
 
320
{
 
321
  int nbPoints = 
 
322
    (serendip ?
 
323
     4 +  6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
 
324
     (order + 1) * (order + 2) * (order + 3) / 6);
 
325
  
 
326
  Double_Matrix point(nbPoints, 3);
 
327
 
 
328
  double overOrder = 1. / order;
 
329
 
 
330
  point(0, 0) = 0.;
 
331
  point(0, 1) = 0.;
 
332
  point(0, 2) = 0.;
 
333
 
 
334
  if (order > 0) {
 
335
    point(1, 0) = order;
 
336
    point(1, 1) = 0;
 
337
    point(1, 2) = 0;
 
338
    
 
339
    point(2, 0) = 0.;
 
340
    point(2, 1) = order;
 
341
    point(2, 2) = 0.;
 
342
    
 
343
    point(3, 0) = 0.;
 
344
    point(3, 1) = 0.;
 
345
    point(3, 2) = order;
 
346
    
 
347
    if (order > 1) {
 
348
      for (int k = 0; k < (order - 1); k++) {
 
349
        point(4 + k, 0) = k + 1;
 
350
        point(4 + order - 1 + k, 0) = order - 1 - k;
 
351
        point(4 + 2 * (order - 1) + k, 0) = 0.; 
 
352
        point(4 + 3 * (order - 1) + k, 0) = 0.;
 
353
        point(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
 
354
        point(4 + 5 * (order - 1) + k, 0) = 0.;
 
355
        
 
356
        point(4 + k, 1) = 0.;
 
357
        point(4 + order - 1 + k, 1) = k + 1;
 
358
        point(4 + 2 * (order - 1) + k, 1) = order - 1 - k; 
 
359
        point(4 + 3 * (order - 1) + k, 1) = 0.;
 
360
        point(4 + 4 * (order - 1) + k, 1) = 0.;
 
361
        point(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
 
362
        
 
363
        point(4 + k, 2) = 0.;
 
364
        point(4 + order - 1 + k, 2) = 0.;
 
365
        point(4 + 2 * (order - 1) + k, 2) = 0.; 
 
366
        point(4 + 3 * (order - 1) + k, 2) = k + 1;
 
367
        point(4 + 4 * (order - 1) + k, 2) = k + 1;
 
368
        point(4 + 5 * (order - 1) + k, 2) = k + 1; 
 
369
      }
 
370
      
 
371
      if (order > 2) {
 
372
        int ns = 4 + 6 * (order - 1);
 
373
        int nbdofface = nbdoftriangle(order - 3);
 
374
        
 
375
        double *u = new double[nbdofface];
 
376
        double *v = new double[nbdofface];
 
377
        double *w = new double[nbdofface];
 
378
        
 
379
        nodepositionface0(order - 3, u, v, w);
 
380
 
 
381
        for (int i = 0; i < nbdofface; i++){
 
382
          point(ns + i, 0) = u[i] * (order - 3) + 1.;
 
383
          point(ns + i, 1) = v[i] * (order - 3) + 1.;
 
384
          point(ns + i, 2) = w[i] * (order - 3);
 
385
        }
 
386
        
 
387
        ns = ns + nbdofface;
 
388
 
 
389
        nodepositionface1(order - 3, u, v, w);
 
390
        
 
391
        for (int i=0; i < nbdofface; i++){
 
392
          point(ns + i, 0) = u[i] * (order - 3) + 1.;
 
393
          point(ns + i, 1) = v[i] * (order - 3) ;
 
394
          point(ns + i, 2) = w[i] * (order - 3) + 1.;
 
395
        }
 
396
        
 
397
        ns = ns + nbdofface;
 
398
 
 
399
        nodepositionface2(order - 3, u, v, w);
 
400
        
 
401
        for (int i = 0; i < nbdofface; i++){
 
402
          point(ns + i, 0) = u[i] * (order - 3) + 1.;
 
403
          point(ns + i, 1) = v[i] * (order - 3) + 1.;
 
404
          point(ns + i, 2) = w[i] * (order - 3) + 1.;
 
405
        }
 
406
 
 
407
        ns = ns + nbdofface;
 
408
 
 
409
        nodepositionface3(order - 3, u, v, w);
 
410
 
 
411
        for (int i = 0; i < nbdofface; i++){
 
412
          point(ns + i, 0) = u[i] * (order - 3);
 
413
          point(ns + i, 1) = v[i] * (order - 3) + 1.;
 
414
          point(ns + i, 2) = w[i] * (order - 3) + 1.;
 
415
        }
 
416
 
 
417
        ns = ns + nbdofface;
 
418
 
 
419
        delete [] u;
 
420
        delete [] v;
 
421
        delete [] w;
 
422
        
 
423
        if (!serendip && order > 3) {
 
424
  
 
425
          Double_Matrix interior = gmshGeneratePointsTetrahedron(order - 4, false);
 
426
          for (int k = 0; k < interior.size1(); k++) {
 
427
            point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
 
428
            point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
 
429
            point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
 
430
          }
 
431
        }
 
432
      }
 
433
    }
 
434
  }
 
435
  
 
436
  point.scale(overOrder);  
 
437
  return point;
 
438
}
 
439
 
 
440
Double_Matrix gmshGeneratePointsTriangle(int order, bool serendip) 
 
441
{
 
442
  int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
 
443
  Double_Matrix point(nbPoints, 2);
 
444
  
 
445
  point(0, 0) = 0;
 
446
  point(0, 1) = 0;
 
447
  
 
448
  double dd = 1. / order;
 
449
 
 
450
  if (order > 0) {
 
451
    point(1, 0) = 1;
 
452
    point(1, 1) = 0;
 
453
    point(2, 0) = 0;
 
454
    point(2, 1) = 1;
 
455
    
 
456
    int index = 3;
 
457
    
 
458
    if (order > 1) {
 
459
      
 
460
      double ksi = 0;
 
461
      double eta = 0;
 
462
      
 
463
      for (int i = 0; i < order - 1; i++, index++) {
 
464
        ksi += dd;
 
465
        point(index, 0) = ksi;
 
466
        point(index, 1) = eta;
 
467
      }
 
468
      
 
469
      ksi = 1.;
 
470
      
 
471
      for (int i = 0; i < order - 1; i++, index++) {
 
472
        ksi -= dd;
 
473
        eta += dd;
 
474
        point(index, 0) = ksi;
 
475
        point(index, 1) = eta;
 
476
      } 
 
477
        
 
478
      eta = 1.;
 
479
      ksi = 0.;
 
480
 
 
481
      for (int i = 0; i < order - 1; i++, index++) {
 
482
        eta -= dd;
 
483
        point(index, 0) = ksi;
 
484
        point(index, 1) = eta;
 
485
      } 
 
486
 
 
487
      if (order > 2 && !serendip) {
 
488
        Double_Matrix inner = gmshGeneratePointsTriangle(order - 3, serendip);
 
489
        inner.scale(1. - 3. * dd);
 
490
        inner.add(dd);
 
491
        Double_Matrix(point.touchSubmatrix(index, nbPoints - index, 0, 2)).memcpy(inner);
 
492
      }
 
493
    }
 
494
  }
 
495
  return point;  
 
496
}
 
497
 
 
498
Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial,
 
499
                                                   const Double_Matrix& point) 
 
500
{
 
501
  if (monomial.size1() != point.size1()) throw;
 
502
  if (monomial.size2() != point.size2()) throw;
 
503
  
 
504
  int ndofs = monomial.size1();
 
505
  int dim   = monomial.size2();
 
506
  
 
507
  Double_Matrix Vandermonde(ndofs, ndofs);
 
508
  
 
509
  for (int i = 0; i < ndofs; i++) {
 
510
    for (int j = 0; j < ndofs; j++) {
 
511
      double dd = 1.;
 
512
      for (int k = 0; k < dim; k++) dd *= pow(point(i, k), monomial(j, k));
 
513
      Vandermonde(i, j) = dd;
 
514
    }
 
515
  }
 
516
  
 
517
  // check for independence
 
518
  
 
519
  double det = Vandermonde.determinant();
 
520
 
 
521
  if (det == 0.0) throw;
 
522
 
 
523
  Double_Matrix coefficient(ndofs, ndofs);
 
524
  
 
525
  for (int i = 0; i < ndofs; i++) {
 
526
    for (int j = 0; j < ndofs; j++) {
 
527
      int f = (i + j) % 2 == 0 ? 1 : -1;
 
528
      Double_Matrix cofactor = Vandermonde.cofactor(i, j);
 
529
      coefficient(i, j) = f * cofactor.determinant() / det;
 
530
    }
 
531
  }
 
532
 
 
533
  Vandermonde.set_all(0.);
 
534
  
 
535
  for (int i = 0; i < ndofs; i++) {
 
536
    for (int j = 0; j < ndofs; j++) {
 
537
      double dd = 1.;
 
538
      for (int k = 0; k < dim; k++) dd *= pow(point(i, k), monomial(j, k));
 
539
      for (int k = 0; k < ndofs; k++) {
 
540
        Vandermonde(i, k) += coefficient(k, j) * dd;
 
541
      }
 
542
    }
 
543
  }
 
544
  return coefficient;
 
545
}
 
546
 
 
547
std::map<int, gmshFunctionSpace> gmshFunctionSpaces::fs;
 
548
 
 
549
const gmshFunctionSpace &gmshFunctionSpaces::find(int tag) 
 
550
{
 
551
  std::map<int,gmshFunctionSpace>::const_iterator it = fs.find(tag);
 
552
  if (it != fs.end()) return it->second;
 
553
  
 
554
  gmshFunctionSpace F;
 
555
  
 
556
  switch (tag){
 
557
  case MSH_TRI_3 :
 
558
    F.monomials = generatePascalTriangle(1);
 
559
    F.points =    gmshGeneratePointsTriangle(1, false);
 
560
    break;
 
561
  case MSH_TRI_6 :
 
562
    F.monomials = generatePascalTriangle(2);
 
563
    F.points =    gmshGeneratePointsTriangle(2, false);
 
564
    break;
 
565
  case MSH_TRI_9 :
 
566
    F.monomials = generatePascalSerendipityTriangle(3);
 
567
    F.points =    gmshGeneratePointsTriangle(3, true);
 
568
    break;
 
569
  case MSH_TRI_10 :
 
570
    F.monomials = generatePascalTriangle(3);
 
571
    F.points =    gmshGeneratePointsTriangle(3, false);
 
572
    break;
 
573
  case MSH_TRI_12 :
 
574
    F.monomials = generatePascalSerendipityTriangle(4);
 
575
    F.points =    gmshGeneratePointsTriangle(4, true);
 
576
    break;
 
577
  case MSH_TRI_15 :
 
578
    F.monomials = generatePascalTriangle(4);
 
579
    F.points =    gmshGeneratePointsTriangle(4, false);
 
580
    break;
 
581
  case MSH_TRI_15I :
 
582
    F.monomials = generatePascalSerendipityTriangle(5);
 
583
    F.points =    gmshGeneratePointsTriangle(5, true);
 
584
    break;
 
585
  case MSH_TRI_21 :
 
586
    F.monomials = generatePascalTriangle(5);
 
587
    F.points =    gmshGeneratePointsTriangle(5, false);
 
588
    break;
 
589
  default :
 
590
    throw;
 
591
  }  
 
592
  F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
 
593
  fs.insert(std::make_pair(tag, F));
 
594
  return fs[tag];
 
595
}