StdProject3D.cpp 18 KB
Newer Older
1 2
#include <cstdio>
#include <cstdlib>
3
#include <cmath>
4

5
#include <boost/format.hpp>
6 7 8
#include <StdRegions/StdHexExp.h>
#include <StdRegions/StdPrismExp.h>
#include <StdRegions/StdNodalPrismExp.h>
9
#include <StdRegions/StdPyrExp.h>
10 11
#include <StdRegions/StdTetExp.h>
#include <StdRegions/StdNodalTetExp.h>
12

13
#include <LibUtilities/Foundations/Foundations.hpp>
14
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
15

16
using namespace Nektar;
17 18 19 20 21
using namespace Nektar::LibUtilities;
using namespace Nektar::StdRegions;

/// Defines a solution which excites all modes in a StdTet expansion.
NekDouble Tet_sol(NekDouble x, NekDouble y, NekDouble z,
22
                  int order1, int order2, int order3);
23

24 25
/// Defines a solution which excites all modes in a StdPrism expansion.
NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z,
26
                    int order1, int order2, int order3);
27

28 29
/// Defines a solution which excites all modes in a StdHex expansion.
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
30 31 32
                  int order1, int order2, int order3,
                  LibUtilities::BasisType btype1, LibUtilities::BasisType btype2,
                  LibUtilities::BasisType btype3);
33 34 35 36

/// This routine projects a polynomial or trigonmetric functions which
/// has energy in all mdoes of the expansions and reports and error
int main(int argc, char *argv[]){
37
    int           i;
38

39 40
    int           order1,order2,order3, nq1,nq2,nq3;
    LibUtilities::PointsType    Qtype1,Qtype2,Qtype3;
41 42 43 44
    LibUtilities::BasisType     btype1 =   LibUtilities::eOrtho_A;
    LibUtilities::BasisType     btype2 =   LibUtilities::eOrtho_B;
    LibUtilities::BasisType     btype3 =   LibUtilities::eOrtho_C;
    LibUtilities::PointsType    NodalType = LibUtilities::eNodalTriElec;
45
    LibUtilities::ShapeType     regionshape;
Dave Moxey's avatar
Dave Moxey committed
46
    StdRegions::StdExpansion *E = NULL;
47 48 49 50 51
    Array<OneD, NekDouble> sol;

    if(argc != 11)
    {
        fprintf(stderr,"Usage: StdProject2D RegionShape Type1 Type2 Type3 "
Dave Moxey's avatar
Dave Moxey committed
52
                       "order1 order2 order3 nq1 nq2 nq3 \n");
53
        fprintf(stderr,"Where RegionShape is an integer value which "
Dave Moxey's avatar
Dave Moxey committed
54
                       "dictates the region shape:\n");
55
        fprintf(stderr,"\t Tetrahedron   = 5\n");
56
        fprintf(stderr,"\t Pyramid       = 6\n");
57 58
        fprintf(stderr,"\t Prism         = 7\n");
        fprintf(stderr,"\t Hexahedron    = 8\n");
59 60

        fprintf(stderr,"Where type is an integer value which "
Dave Moxey's avatar
Dave Moxey committed
61
                       "dictates the basis as:\n");
62

63 64 65 66 67 68 69 70
        fprintf(stderr,"\t Ortho_A             =  1\n");
        fprintf(stderr,"\t Ortho_B             =  2\n");
        fprintf(stderr,"\t Ortho_C             =  3\n");
        fprintf(stderr,"\t Modified_A          =  4\n");
        fprintf(stderr,"\t Modified_B          =  5\n");
        fprintf(stderr,"\t Modified_C          =  6\n");
        fprintf(stderr,"\t Fourier             =  7\n");
        fprintf(stderr,"\t Lagrange            =  8\n");
Daniele de Grazia's avatar
Daniele de Grazia committed
71
        fprintf(stderr,"\t Gauss Lagrange      =  9\n");
72 73 74 75 76 77 78
        fprintf(stderr,"\t Legendre            = 10\n");
        fprintf(stderr,"\t Chebyshev           = 11\n");
        fprintf(stderr,"\t Nodal tri (Electro) = 12\n");
        fprintf(stderr,"\t Nodal tri (Fekete)  = 13\n");
        fprintf(stderr,"\t Nodal tet (Electro) = 14\n");
        fprintf(stderr,"\t Nodal tet (Even)    = 15\n");
        fprintf(stderr,"\t Nodal prism (Even)  = 16\n");
79 80 81 82

        exit(1);
    }

83
    regionshape = (LibUtilities::ShapeType) atoi(argv[1]);
84 85

    // Check to see if 3D region
86
    if ((regionshape != LibUtilities::eTetrahedron)
87
        && (regionshape != LibUtilities::ePyramid)
88 89
        && (regionshape != LibUtilities::ePrism)
        && (regionshape != LibUtilities::eHexahedron))
90 91 92 93 94 95 96
    {
        NEKERROR(ErrorUtil::efatal,"This shape is not a 3D region");
    }

    int btype1_val = atoi(argv[2]);
    int btype2_val = atoi(argv[3]);
    int btype3_val = atoi(argv[4]);
97
    
98
    if (btype1_val <= 11 && btype2_val <= 11)
99 100 101 102 103
    {
        btype1 =   (LibUtilities::BasisType) btype1_val;
        btype2 =   (LibUtilities::BasisType) btype2_val;
        btype3 =   (LibUtilities::BasisType) btype3_val;
    }
104
    else if(btype1_val >=12 && btype2_val <= 16)
105
    {
106
        if (regionshape == LibUtilities::eTetrahedron)
107 108 109 110 111
        {
            btype1 = LibUtilities::eOrtho_A;
            btype2 = LibUtilities::eOrtho_B;
            btype3 = LibUtilities::eOrtho_C;
        }
112
        else if (regionshape == LibUtilities::ePrism)
113 114 115 116 117 118
        {
            btype1 = LibUtilities::eOrtho_A;
            btype2 = LibUtilities::eOrtho_A;
            btype3 = LibUtilities::eOrtho_B;
        }
        
119
        if(btype1_val == 12)
120 121 122
        {
            NodalType = LibUtilities::eNodalTriElec;
        }
123
        else if (btype1_val == 13)
124 125 126
        {
            NodalType = LibUtilities::eNodalTriFekete;
        }
127
        else if (btype1_val == 14)
128 129 130
        {
            NodalType = LibUtilities::eNodalTetElec;
        }
131
        else if (btype1_val == 15)
132 133 134 135 136 137 138
        {
            NodalType = LibUtilities::eNodalTetEvenlySpaced;
        }
        else
        {
            NodalType = LibUtilities::eNodalPrismEvenlySpaced;
        }
139 140 141 142 143
    }

    // Check to see that correct Expansions are used
    switch(regionshape)
    {
144
        case LibUtilities::eTetrahedron:
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
            if((btype1 == eOrtho_B) || (btype1 == eOrtho_C)
               || (btype1 == eModified_B) || (btype1 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 1 cannot be of type Ortho_B, "
                         "Ortho_C, Modified_B or Modified_C");
            }
            if((btype2 == eOrtho_A) || (btype2 == eOrtho_C)
               || (btype2 == eModified_A) || (btype2 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 2 cannot be of type Ortho_A, "
                         "Ortho_C, Modified_A or Modified_C");
            }
            if((btype3 == eOrtho_A) || (btype3 == eOrtho_B)
               || (btype3 == eModified_A) || (btype3 == eModified_B))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 3 cannot be of type Ortho_A, "
                         "Ortho_B, Modified_A or Modified_B");
            }
            break;
164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
        case LibUtilities::ePyramid:
            if((btype1 == eOrtho_B) || (btype1 == eOrtho_C)
               || (btype1 == eModified_B) || (btype1 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 1 cannot be of type Ortho_B, "
                         "Ortho_C, Modified_B or Modified_C");
            }
            if((btype2 == eOrtho_B) || (btype2 == eOrtho_C)
               || (btype2 == eModified_B) || (btype2 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 2 cannot be of type Ortho_B, "
                         "Ortho_C, Modified_B or Modified_C");
            }
            if((btype3 == eOrtho_A) || (btype3 == eOrtho_B)
               || (btype3 == eModified_A) || (btype3 == eModified_B))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 3 cannot be of type Ortho_A, "
                         "Ortho_B, Modified_A or Modified_B");
            }
            break;
184
        case LibUtilities::ePrism:
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
            if((btype1 == eOrtho_B) || (btype1 == eOrtho_C)
               || (btype1 == eModified_B) || (btype1 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 1 cannot be of type Ortho_B, "
                         "Ortho_C, Modified_B or Modified_C");
            }
            if((btype2 == eOrtho_B) || (btype2 == eOrtho_C)
               || (btype2 == eModified_B) || (btype2 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 2 cannot be of type Ortho_B, "
                         "Ortho_C, Modified_B or Modified_C");
            }
            if((btype3 == eOrtho_A) || (btype3 == eOrtho_C)
               || (btype3 == eModified_A) || (btype3 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 3 cannot be of type Ortho_A, "
                         "Ortho_C, Modified_A or Modified_C");
            }
            break;
204
        case LibUtilities::eHexahedron:
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
            if((btype1 == eOrtho_B) || (btype1 == eOrtho_C)
               || (btype1 == eModified_B) || (btype1 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 1 is for 2 or 3D expansions");
            }
            if((btype2 == eOrtho_B) || (btype2 == eOrtho_C)
               || (btype2 == eModified_B) || (btype2 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 2 is for 2 or 3D expansions");
            }
            if((btype3 == eOrtho_B) || (btype3 == eOrtho_C)
               || (btype3 == eModified_B) || (btype3 == eModified_C))
            {
                NEKERROR(ErrorUtil::efatal, "Basis 3 is for 2 or 3D expansions");
            }
            break;
        default:
            ASSERTL0(false, "Not a 3D expansion.");
            break;
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
    }

    order1 =   atoi(argv[5]);
    order2 =   atoi(argv[6]);
    order3 =   atoi(argv[7]);
    nq1    =   atoi(argv[8]);
    nq2    =   atoi(argv[9]);
    nq3    =   atoi(argv[10]);

    sol = Array<OneD, NekDouble>(nq1*nq2*nq3);

    if(btype1 != LibUtilities::eFourier)
    {
        Qtype1 = LibUtilities::eGaussLobattoLegendre;
    }
    else
    {
        Qtype1 = LibUtilities::eFourierEvenlySpaced;
    }

    if(btype2 != LibUtilities::eFourier)
    {
246
        if (regionshape == LibUtilities::eTetrahedron) 
247
        {
248 249 250 251 252 253
            Qtype2 = LibUtilities::eGaussRadauMAlpha1Beta0;
        }
        else
        {
            Qtype2 = LibUtilities::eGaussLobattoLegendre;
        }
254 255 256 257 258 259 260 261
    }
    else
    {
        Qtype2 = LibUtilities::eFourierEvenlySpaced;
    }

    if(btype3 != LibUtilities::eFourier)
    {
262 263
        if (regionshape == LibUtilities::eTetrahedron ||
            regionshape == LibUtilities::ePyramid) 
264
        {
265 266
            Qtype3 = LibUtilities::eGaussRadauMAlpha2Beta0;
        }
267
        else if (regionshape == LibUtilities::ePrism)
268 269 270
        {
            Qtype3 = LibUtilities::eGaussRadauMAlpha1Beta0;
        }
271 272 273 274
        else
        {
            Qtype3 = LibUtilities::eGaussLobattoLegendre;
        }
275 276 277 278 279 280 281 282 283 284 285 286 287 288
    }
    else
    {
        Qtype3 = LibUtilities::eFourierEvenlySpaced;
    }

    //-----------------------------------------------
    // Define a 3D expansion based on basis definition
    Array<OneD,NekDouble> x(nq1*nq2*nq3);
    Array<OneD,NekDouble> y(nq1*nq2*nq3);
    Array<OneD,NekDouble> z(nq1*nq2*nq3);

    switch(regionshape)
    {
289
    case LibUtilities::eTetrahedron:
290 291 292 293 294 295 296 297
        {
            const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
            const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
            const LibUtilities::PointsKey Pkey3(nq3,Qtype3);
            const LibUtilities::BasisKey  Bkey1(btype1,order1,Pkey1);
            const LibUtilities::BasisKey  Bkey2(btype2,order2,Pkey2);
            const LibUtilities::BasisKey  Bkey3(btype3,order3,Pkey3);

298 299 300 301 302 303 304 305
            if(btype1_val >= 10)
            {
                E = new StdRegions::StdNodalTetExp(Bkey1,Bkey2,Bkey3,NodalType);
            }
            else
            {
                E = new StdRegions::StdTetExp(Bkey1,Bkey2,Bkey3);
            }
306

307
            E->GetCoords(x,y,z);
308
            
309 310 311 312 313 314 315 316 317
            //----------------------------------------------
            // Define solution to be projected
            for(i = 0; i < nq1*nq2*nq3; ++i)
            {
                sol[i]  = Tet_sol(x[i],y[i],z[i],order1,order2,order3);
            }
            //----------------------------------------------
        }
        break;
318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334
    case LibUtilities::ePyramid:
        {
            const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
            const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
            const LibUtilities::PointsKey Pkey3(nq3,Qtype3);
            const LibUtilities::BasisKey  Bkey1(btype1,order1,Pkey1);
            const LibUtilities::BasisKey  Bkey2(btype2,order2,Pkey2);
            const LibUtilities::BasisKey  Bkey3(btype3,order3,Pkey3);

            StdRegions::StdPyrExp *F = new StdRegions::StdPyrExp(Bkey1,Bkey2,Bkey3);
            E = F;
            E->GetCoords(x,y,z);

            //----------------------------------------------
            // Define solution to be projected
            for(i = 0; i < nq1*nq2*nq3; ++i)
            {
Dave Moxey's avatar
Dave Moxey committed
335
                sol[i]  = Tet_sol(x[i],y[i],z[i],order1,order2,order3);
336 337 338 339
            }
            //----------------------------------------------
        }
        break;
340
    case LibUtilities::ePrism:
341 342 343 344 345 346 347 348
        {
            const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
            const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
            const LibUtilities::PointsKey Pkey3(nq3,Qtype3);
            const LibUtilities::BasisKey  Bkey1(btype1,order1,Pkey1);
            const LibUtilities::BasisKey  Bkey2(btype2,order2,Pkey2);
            const LibUtilities::BasisKey  Bkey3(btype3,order3,Pkey3);

349 350 351 352 353 354 355 356
            if (btype1_val >= 10)
            {
                E = new StdRegions::StdNodalPrismExp(Bkey1,Bkey2,Bkey3,NodalType);
            }
            else
            {
                E = new StdRegions::StdPrismExp(Bkey1,Bkey2,Bkey3);
            }
357 358 359 360 361 362 363 364 365 366 367
            E->GetCoords(x,y,z);

            //----------------------------------------------
            // Define solution to be projected
            for(i = 0; i < nq1*nq2*nq3; ++i)
            {
                sol[i]  = Prism_sol(x[i],y[i],z[i],order1,order2,order3);
            }
            //----------------------------------------------
        }
        break;
368
    case LibUtilities::eHexahedron:
369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
        {
            const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
            const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
            const LibUtilities::PointsKey Pkey3(nq3,Qtype3);
            const LibUtilities::BasisKey Bkey1(btype1,order1,Pkey1);
            const LibUtilities::BasisKey Bkey2(btype2,order2,Pkey2);
            const LibUtilities::BasisKey Bkey3(btype3,order3,Pkey3);
            E = new StdRegions::StdHexExp(Bkey1,Bkey2,Bkey3);

            //----------------------------------------------
            // Define solution to be projected
            Array<OneD,NekDouble> x = Array<OneD,NekDouble>(nq1*nq2*nq3);
            Array<OneD,NekDouble> y = Array<OneD,NekDouble>(nq1*nq2*nq3);
            Array<OneD,NekDouble> z = Array<OneD,NekDouble>(nq1*nq2*nq3);
            E->GetCoords(x,y,z);

            for(i = 0; i < nq1*nq2*nq3; ++i)
            {
                sol[i]  = Hex_sol(x[i],y[i],z[i],order1,order2,order3,
                                  btype1,btype2,btype3);
            }
            //---------------------------------------------
        }
        break;
393 394 395
        default:
            ASSERTL0(false, "Not a 3D expansion.");
            break;
396 397
    }

398 399 400 401
    
    Array<OneD, NekDouble> phys (nq1*nq2*nq3);
    Array<OneD, NekDouble> coeffs(order1*order2*order3);

402 403
    //---------------------------------------------
    // Project onto Expansion
404
    E->FwdTrans(sol,coeffs);
405 406 407 408
    //---------------------------------------------

    //-------------------------------------------
    // Backward Transform Solution to get projected values
409
    E->BwdTrans(coeffs,phys);
410 411 412 413
    //-------------------------------------------

    //--------------------------------------------
    // Calculate L_inf error
414 415
    cout << "L infinity error: " << E->Linf(phys,sol) << endl;
    cout << "L 2 error:        " << E->L2  (phys,sol) << endl;
416 417 418 419 420 421 422 423 424
    //--------------------------------------------

    //-------------------------------------------
    // Evaulate solution at x = y = z = 0  and print error
    Array<OneD, NekDouble> t = Array<OneD, NekDouble>(3);
    t[0] = -0.5;
    t[1] = -0.25;
    t[2] = -0.3;

425 426
    if(regionshape == LibUtilities::eTetrahedron ||
       regionshape == LibUtilities::ePyramid)
427
    {
Dave Moxey's avatar
Dave Moxey committed
428
        sol[0] = Tet_sol(t[0], t[1], t[2], order1, order2, order3);
429
    }
430
    else if (regionshape == LibUtilities::ePrism)
431 432 433
    {
        sol[0] = Prism_sol(t[0], t[1], t[2], order1, order2, order3);
    }
434 435 436 437 438 439
    else
    {
        sol[0] = Hex_sol(t[0], t[1], t[2], order1, order2, order3,
                         btype1, btype2, btype3);
    }

440
    NekDouble nsol = E->PhysEvaluate(x,phys);
441 442 443 444 445 446 447 448 449
    cout << "error at x = (" << t[0] << "," << t[1] << "," << t[2] << "): "
         << nsol - sol[0] << endl;
    //-------------------------------------------

    return 0;
}


NekDouble Tet_sol(NekDouble x, NekDouble y, NekDouble z,
450
                  int order1, int order2, int order3)
451 452
{
    int    l,k,m;
453
    NekDouble sol = 0.0;
454 455 456 457 458 459 460

    for(k = 0; k < order1; ++k)
    {
        for(l = 0; l < order2-k; ++l)
        {
            for(m = 0; m < order3-k-l; ++m)
            {
Dave Moxey's avatar
Dave Moxey committed
461
                sol += pow(x,k)*pow(y,l)*pow(z,m);
462 463 464
            }
        }
    }
Dave Moxey's avatar
Dave Moxey committed
465

466 467 468
    return sol;
}

469 470 471 472 473 474 475 476 477 478 479 480
NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z,
                    int order1, int order2, int order3)
{
    int    l,k,m;
    NekDouble sol = 0;

    for(k = 0; k < order1; ++k)
    {
        for(l = 0; l < order2; ++l)
        {
            for(m = 0; m < order3-k; ++m)
            {
481
                sol += pow(x,k)*pow(y,l)*pow(z,m);;
482 483 484 485 486 487
            }
        }
    }
    return sol;
}

488
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
489 490 491 492
                  int order1, int order2, int order3,
                  LibUtilities::BasisType btype1,
                  LibUtilities::BasisType btype2,
                  LibUtilities::BasisType btype3)
493
{
494
    int i,j,k;
495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520
    NekDouble sol = 0.0;

    int  Nx = (btype1 == LibUtilities::eFourier ? order1/2 : order1);
    int  Ny = (btype2 == LibUtilities::eFourier ? order2/2 : order2);
    int  Nz = (btype3 == LibUtilities::eFourier ? order3/2 : order3);
    bool Fx = (btype1 == LibUtilities::eFourier);
    bool Fy = (btype2 == LibUtilities::eFourier);
    bool Fz = (btype3 == LibUtilities::eFourier);
    NekDouble a;

    for (i = 0; i < Nx; ++i)
    {
        for (j = 0; j < Ny; ++j)
        {
            for (k = 0; k < Nz; ++k)
            {
                a  = (Fx ? sin(M_PI*i*x) + cos(M_PI*i*x) : pow(x,i));
                a *= (Fy ? sin(M_PI*j*y) + cos(M_PI*j*y) : pow(y,j));
                a *= (Fz ? sin(M_PI*k*z) + cos(M_PI*k*z) : pow(z,k));
                sol += a;
            }
        }
    }
    return sol;
}