StdProject_Diff2D.cpp 13.1 KB
Newer Older
1
2
#include <cstdio>
#include <cstdlib>
3
4
5
#include <StdRegions/StdQuadExp.h>
#include <StdRegions/StdTriExp.h>
#include <StdRegions/StdNodalTriExp.h>
6

7
#include <LibUtilities/Foundations/Foundations.hpp>
8
9
10
11
12

using namespace Nektar;

NekDouble Tri_sol(NekDouble x, NekDouble y, int order1, int order2);
NekDouble Tri_Dsol(NekDouble x, NekDouble y, int order1, int order2);
13
NekDouble Quad_sol(NekDouble x, NekDouble y, int order1, int order2,
14
                   LibUtilities::BasisType btype1, LibUtilities::BasisType btype2);
15
NekDouble Quad_Dsol(NekDouble x, NekDouble y, int order1, int order2,
16
                    LibUtilities::BasisType btype1, LibUtilities::BasisType btype2);
17

18
// This routine projects a polynomial or trigonmetric functions which
19
20
21
22
23
// has energy in all mdoes of the expansions and reports and error
//

static double  pow_loc(const double val, const int i)
{
24
    return (i < 0)? 1.0: pow(val,i);
25
26
}

Blake Nelson's avatar
Blake Nelson committed
27
int main(int argc, char *argv[])
28
{
29
30
31
32
    int           i;

    int           order1,order2, nq1,nq2;
    LibUtilities::PointsType    Qtype1,Qtype2;
33
34
35
    LibUtilities::BasisType     btype1 =   LibUtilities::eOrtho_A;
    LibUtilities::BasisType     btype2 =   LibUtilities::eOrtho_B;
    LibUtilities::PointsType    NodalType = LibUtilities::eNodalTriElec;
36
    LibUtilities::ShapeType     regionshape;
37
38
39
40
41
42
43
44
45
46
    StdRegions::StdExpansion          *E;
    Array<OneD, NekDouble>  sol,dx,dy,x,y;

    if(argc != 8)
    {
        fprintf(stderr,"Usage: StdProject2D_Diff RegionShape Type1 Type2 order1 "
                "order2  nq1 nq2  \n");

        fprintf(stderr,"Where RegionShape is an integer value which "
                "dictates the region shape:\n");
47
48
        fprintf(stderr,"\t Triangle      = 3\n");
        fprintf(stderr,"\t Quadrilateral = 4\n");
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69


        fprintf(stderr,"Where type is an integer value which "
                "dictates the basis as:\n");

        fprintf(stderr,"\t Ortho_A    = 1\n");
        fprintf(stderr,"\t Ortho_B    = 2\n");
        fprintf(stderr,"\t Modified_A = 4\n");
        fprintf(stderr,"\t Modified_B = 5\n");
        fprintf(stderr,"\t Fourier    = 7\n");
        fprintf(stderr,"\t Lagrange   = 8\n");
        fprintf(stderr,"\t Legendre   = 9\n");
        fprintf(stderr,"\t Chebyshev  = 10\n");
        fprintf(stderr,"\t Nodal tri (Electro) = 11\n");
        fprintf(stderr,"\t Nodal tri (Fekete)  = 12\n");

        fprintf(stderr,"Note type = 3,6 are for three-dimensional basis\n");

        exit(1);
    }

70
    regionshape = (LibUtilities::ShapeType) atoi(argv[1]);
71
72

    // Check to see if 2D region
73
    if((regionshape != LibUtilities::eTriangle)&&(regionshape != LibUtilities::eQuadrilateral))
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
    {
        NEKERROR(ErrorUtil::efatal,"This shape is not a 2D region");
    }

    int btype1_val =  atoi(argv[2]);
    int btype2_val =  atoi(argv[3]);

    if(( btype1_val <= 10)&&( btype2_val <= 10))
    {
        btype1 =   (LibUtilities::BasisType) btype1_val;
        btype2 =   (LibUtilities::BasisType) btype2_val;
    }
    else if(( btype1_val >=11)&&(btype2_val <= 12))
    {
        btype1 =   LibUtilities::eOrtho_A;
        btype2 =   LibUtilities::eOrtho_B;

        if(btype1_val == 11)
        {
            NodalType = LibUtilities::eNodalTriElec;
        }
        else
        {
            NodalType = LibUtilities::eNodalTriFekete;
        }
    }


    // Check to see that correct Expansions are used
    switch(regionshape)
104
    {
105
106
107
108
109
110
    case LibUtilities::eTriangle:
        if((btype1 == LibUtilities::eOrtho_B)||(btype1 == LibUtilities::eModified_B))
        {
            NEKERROR(ErrorUtil::efatal,
                     "Basis 1 cannot be of type Ortho_B or Modified_B");
        }
111

112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
        if((btype2 != LibUtilities::eOrtho_B)&&(btype2 != LibUtilities::eModified_B))
        {
            NEKERROR(ErrorUtil::efatal,
                     "Basis 2 must be of type Ortho_B or Modified_B");
        }
        break;
    case LibUtilities::eQuadrilateral:
        if((btype1 == LibUtilities::eOrtho_B)||(btype1 == LibUtilities::eOrtho_B)||
           (btype1 == LibUtilities::eModified_B)||(btype1 == LibUtilities::eModified_C))
        {
            NEKERROR(ErrorUtil::efatal,
                     "Basis 1 is for 2 or 3D expansions");
        }
        
        if((btype2 == LibUtilities::eOrtho_B)||(btype2 == LibUtilities::eOrtho_B)||
           (btype2 == LibUtilities::eModified_B)||(btype2 == LibUtilities::eModified_C))
        {
            NEKERROR(ErrorUtil::efatal, "Basis 2 is for 2 or 3D expansions");
        }
        break;
        
    default:
        ASSERTL0(false, "Not a 2D expansion.");
        break;
136
    }
137
    
138
139
140
141
142
143
144
145
146
147
148
149
    order1 =   atoi(argv[4]);
    order2 =   atoi(argv[5]);
    nq1    =   atoi(argv[6]);
    nq2    =   atoi(argv[7]);

    dx  = Array<OneD, NekDouble>(nq1*nq2);
    dy  = Array<OneD, NekDouble>(nq1*nq2);
    sol = Array<OneD, NekDouble>(nq1*nq2);
    x   = Array<OneD, NekDouble>(nq1*nq2);
    y   = Array<OneD, NekDouble>(nq1*nq2);

    if(btype1 != LibUtilities::eFourier)
150
    {
151
        Qtype1 = LibUtilities::eGaussLobattoLegendre;
152
    }
153
154
155
156
157
158
159
    else
    {
        Qtype1 = LibUtilities::eFourierEvenlySpaced;
    }

    if(btype2 != LibUtilities::eFourier)
    {
160
161
        if (regionshape == LibUtilities::eTriangle) 
        {
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
            Qtype2 = LibUtilities::eGaussRadauMAlpha1Beta0;
        }
        else
        {
            Qtype2 = LibUtilities::eGaussLobattoLegendre;
        }
    }
    else
    {
        Qtype2 = LibUtilities::eFourierEvenlySpaced;
    }


    //-----------------------------------------------
    // Define a segment expansion based on basis definition

    switch(regionshape)
    {
180
        case LibUtilities::eTriangle:
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
        {
            const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
            const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
            const LibUtilities::BasisKey  Bkey1(btype1,order1,Pkey1);
            const LibUtilities::BasisKey  Bkey2(btype2,order2,Pkey2);


            if(btype1_val >= 10)
            {
                E = new StdRegions::StdNodalTriExp(Bkey1,Bkey2,NodalType);
            }
            else
            {
                E = new StdRegions::StdTriExp(Bkey1,Bkey2);
            }

            E->GetCoords(x,y);

            //----------------------------------------------
            // Define solution to be differentiated
            for(i = 0; i < nq1*nq2; ++i)
            {
                sol[i] = Tri_sol(x[i],y[i],order1,order2);
            }
            //----------------------------------------------
        }
        break;
208
        case LibUtilities::eQuadrilateral:
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
        {
            const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
            const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
            const LibUtilities::BasisKey  Bkey1(btype1,order1,Pkey1);
            const LibUtilities::BasisKey  Bkey2(btype2,order2,Pkey2);
            E = new StdRegions::StdQuadExp (Bkey1,Bkey2);

            //----------------------------------------------
            // Define solution to be differentiated
            E->GetCoords(x,y);

            for(i = 0; i < nq1*nq2; ++i)
            {
                sol[i]  = Quad_sol(x[i],y[i],order1,order2,btype1,btype2);
            }
            //---------------------------------------------
        }
        break;
        default:
            ASSERTL0(false, "Not a 2D expansion.");
            break;
    }

    //---------------------------------------------
    // Evaluate derivative of solution, add together and put in sol
    E->PhysDeriv(sol,dx,dy);
    Vmath::Vadd(nq1*nq2,dx,1,dy,1,sol,1);
    //---------------------------------------------

    //---------------------------------------------
    // Project onto Expansion
    E->FwdTrans(sol,E->UpdateCoeffs());
    //---------------------------------------------

    //-------------------------------------------
    // Backward Transform Solution to get projected values
    E->BwdTrans(E->GetCoeffs(),E->UpdatePhys());
    //-------------------------------------------

    //----------------------------------------------
    // Define exact solution of differential
    switch(regionshape)
    {
252
        case LibUtilities::eTriangle:
253
254
255
256
257
258
259
260
261
262
        {
            //----------------------------------------------
            // Define solution to be differentiated
            for(i = 0; i < nq1*nq2; ++i)
            {
                sol[i] = Tri_Dsol(x[i],y[i],order1,order2);
            }
            //----------------------------------------------
        }
        break;
263
        case LibUtilities::eQuadrilateral:
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
        {
            for(i = 0; i < nq1*nq2; ++i)
            {
                sol[i] = Quad_Dsol(x[i],y[i],order1,order2,btype1,btype2);
            }
        }
        //---------------------------------------------
        break;
        default:
            ASSERTL0(false, "Not a 2D expansion.");
            break;
    }

    //--------------------------------------------
    // Calculate L_inf error
    cout << "L infinity error: " << E->Linf(sol) << endl;
    cout << "L 2 error:        " << E->L2  (sol) << endl;
    //--------------------------------------------
    return 0;
283
284
}

285
286
NekDouble Tri_sol(NekDouble x, NekDouble y, int order1, int order2)
{
287
288
    int    l,k;
    NekDouble sol = 0;
289

290
291
292
293
294
295
296
    for(k = 0; k < order1; ++k)
    {
        for(l = 0; l < order2-k; ++l)
        {
            sol += pow(x,k)*pow(y,l);
        }
    }
297

298
299
300
    return sol;
}

301
302
NekDouble Tri_Dsol(NekDouble x, NekDouble y, int order1, int order2)
{
303
304
    int    l,k;
    NekDouble sol = 0;
305

306
307
308
309
    for(k = 0; k < order1; ++k)
    {
        for(l = 0; l < order2-k; ++l)
        {
310
            sol +=  k*pow_loc(x,k-1)*pow_loc(y,l) +
311
312
313
                l*pow_loc(x,k)*pow_loc(y,l-1);
        }
    }
314

315
316
317
    return sol;
}

318
319
NekDouble Quad_sol(NekDouble x, NekDouble y, int order1, int order2, LibUtilities::BasisType btype1, LibUtilities::BasisType btype2)
{
320

321
322
    int k,l;
    NekDouble sol = 0;
323

324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    if(btype1 != LibUtilities::eFourier)
    {
        if(btype2 != LibUtilities::eFourier)
        {
            for(k = 0; k < order1; ++k)
            {
                for(l = 0; l < order2; ++l)
                {
                    sol += pow(x,k)*pow(y,l);
                }
            }
        }
        else
        {
            for(k = 0; k < order1; ++k)
            {
                for(l = 0; l < order2/2; ++l)
                {
                    sol += pow(x,k)*sin(M_PI*l*y) + pow(x,k)*cos(M_PI*l*y);
                }
            }
        }
    }
    else
    {
        if(btype2 != LibUtilities::eFourier)
        {
            for(k = 0; k < order1/2; ++k)
            {
                for(l = 0; l < order2; ++l)
                {
                    sol += sin(M_PI*k*x)*pow(y,l) + cos(M_PI*k*x)*pow(y,l);
                }
            }
        }
        else
        {
            for(k = 0; k < order1/2; ++k)
            {
                for(l = 0; l < order2/2; ++l)
                {
                    sol += sin(M_PI*k*x)*sin(M_PI*l*y)
                        + sin(M_PI*k*x)*cos(M_PI*l*y)
                        + cos(M_PI*k*x)*sin(M_PI*l*y)
                        + cos(M_PI*k*x)*cos(M_PI*l*y);
                }
            }
        }
    }
373

374
375
376
377
    return sol;
}


378
379
NekDouble Quad_Dsol(NekDouble x, NekDouble y, int order1, int order2, LibUtilities::BasisType btype1,LibUtilities::BasisType btype2)
{
380

381
382
    int k,l;
    NekDouble sol = 0;
383

384
385
386
387
388
389
390
391
    if(btype1 != LibUtilities::eFourier)
    {
        if(btype2 !=LibUtilities:: eFourier)
        {
            for(k = 0; k < order1; ++k)
            {
                for(l = 0; l < order2; ++l)
                {
392
                    sol +=  k*pow_loc(x,k-1)*pow_loc(y,l)
393
394
395
396
397
398
399
400
401
402
                        + l*pow_loc(x,k)*pow_loc(y,l-1);
                }
            }
        }
        else
        {
            for(k = 0; k < order1; ++k)
            {
                for(l = 0; l < order2/2; ++l)
                {
403
404
405
                    sol += k*pow_loc(x,k-1)*sin(M_PI*l*y)
                        + M_PI*l*pow_loc(x,k)*cos(M_PI*l*y) +
                        + k*pow_loc(x,k-1)*cos(M_PI*l*y)
406
407
408
409
410
411
412
413
414
415
416
417
418
                        - M_PI*l*pow_loc(x,k)*sin(M_PI*l*y);
                }
            }
        }
    }
    else
    {
        if(btype2 != LibUtilities::eFourier)
        {
            for(k = 0; k < order1/2; ++k)
            {
                for(l = 0; l < order2; ++l)
                {
419
420
421
                    sol += M_PI*k*cos(M_PI*k*x)*pow_loc(y,l)
                        + l*sin(M_PI*k*x)*pow_loc(y,l-1) +
                        - M_PI*k*sin(M_PI*k*x)*pow_loc(y,l)
422
423
424
425
426
427
428
429
430
431
432
                        + l*sin(M_PI*k*x)*pow_loc(y,l-1);
                }
            }
        }
        else
        {
            for(k = 0; k < order1/2; ++k)
            {
                for(l = 0; l < order2/2; ++l)
                {
                    sol += M_PI*k*cos(M_PI*k*x)*sin(M_PI*l*y)
433
434
435
436
437
438
439
                        +  M_PI*l*sin(M_PI*k*x)*cos(M_PI*l*y)
                        +  M_PI*k*cos(M_PI*k*x)*cos(M_PI*l*y)
                        -  M_PI*l*sin(M_PI*k*x)*sin(M_PI*l*y)
                        -  M_PI*k*sin(M_PI*k*x)*sin(M_PI*l*y)
                        +  M_PI*l*cos(M_PI*k*x)*cos(M_PI*l*y)
                        -  M_PI*k*sin(M_PI*k*x)*cos(M_PI*l*y)
                        -  M_PI*l*cos(M_PI*k*x)*sin(M_PI*l*y);
440
441
442
443
                }
            }
        }
    }
444

445
    return sol;
446
}