Cartesian.cpp 13.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
#include "gtest/gtest.h"

#include "opal_test_utilities/SilenceTest.h"

#include "Meshes/Cartesian.h"

#include "FieldLayout/CenteredFieldLayout.h"

constexpr unsigned int DIM = 3;
constexpr double roundOffError = 1e-10;

TEST(Meshes, Cartesian)
{
    OpalTestUtilities::SilenceTest silencer;

    const unsigned D = DIM; // Hardwire dimensionality
    const unsigned nv = 6;  // Hardwire number of vertices in every direction
    unsigned vnodes = 4;    // Hardwire 4 vnodes

    // Sizes:
    unsigned nverts[D], ncells[D];
22
    unsigned totverts = 1, totcells = 1;
23 24
    unsigned int d;

25
    for (d = 0; d < D; d++) {
26 27 28 29 30 31
        ncells[d] = nv - 1;
        nverts[d] = nv;
        totcells *= ncells[d];
        totverts *= nverts[d];
    }
    NDIndex<D> verts, cells;
32
    for (d = 0; d < D; d++) {
33 34 35 36 37 38 39 40 41 42 43 44
        verts[d] = Index(nverts[d]);
        cells[d] = Index(ncells[d]);
    }

    //---------------------------------------------------------------------------
    // Construct some CenteredFieldLayout's and Field's to be used below:

    // Create cartesian mesh object:
    typedef Cartesian<D,double> M;

    double* delX[D];

45
    for (d = 0; d < D; d++)
46 47 48
        delX[d] = new double[nverts[d]];

    Vektor<double,D> origin;
49
    for (d = 0; d < D; d++)
50 51 52
        origin(d) = d + 1.0;

    // Assign nonuniform mesh-spacing values to each component (linear ramps):
53
    for (d = 0; d < D; d++) {
54
        double multipplier = (d + 1)*1.0;
55
        for (unsigned int vert = 0; vert < nverts[d]; vert++) {
56 57 58 59 60 61
            (delX[d])[vert] = multipplier*(1 + vert);
        }
    }

    // Mesh boundary conditions:
    MeshBC_E mbc[2*D];
62
    for (unsigned b = 0; b < (2*D); b++)
63 64 65 66 67 68 69 70 71
        mbc[b] = Reflective;

    // Test constructing mesh, and then setting spacing, origin, BC's
    M mesh(verts);
    mesh.set_origin(origin);
    mesh.set_meshSpacing(delX);
    mesh.set_MeshBC(mbc);

    // Clean up mesh spacing arrays
72
    for (d = 0; d < D; d++)
73 74 75 76 77 78 79
        delete [] delX[d];

    // ada have to cross check Div() fails without this
    mesh.storeSpacingFields();

    // Construct CenteredFieldLayout's using this for Vert and Cell centering:
    e_dim_tag edt[D];
80
    for (d = 0; d < D; d++) edt[d] = PARALLEL;
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
    CenteredFieldLayout<D,M,Cell> cl(mesh, edt, vnodes);
    CenteredFieldLayout<D,M,Vert> vl(mesh, edt, vnodes);

    // Use 1 guard layer in all Field's:
    GuardCellSizes<D> gc(1);

    // Vectors:
    BConds<Vektor<double,D>,D,M,Vert> vvbc;
    BConds<Vektor<double,D>,D,M,Cell> vcbc;

    // Scalars:
    BConds<double,D,M,Cell> scbc;

    // Symmetric tensors:
    BConds<SymTenzor<double,D>,D,M,Cell> stcbc;

    // Tensors:
    BConds<Tenzor<double,D>,D,M,Cell> tcbc;

    // Use linear negative reflecting conditions:
101
    for (unsigned int face = 0; face < 2*D; face++) {
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
        vvbc[face]  = new NegReflectFace<Vektor<double,D>,D,M,Vert>(face);
        vcbc[face]  = new NegReflectFace<Vektor<double,D>,D,M,Cell>(face);
        scbc[face]  = new NegReflectFace<double,D,M,Cell>(face);
        stcbc[face] = new NegReflectFace<SymTenzor<double,D>,D,M,Cell>(face);
        tcbc[face] =  new NegReflectFace<Tenzor<double,D>,D,M,Cell>(face);
    }

    // Now use all this to construct some Field's:
    Field<Vektor<double,D>,D,M,Vert> vectorVert(mesh, vl, gc, vvbc);
    Field<Vektor<double,D>,D,M,Cell> vectorCell(mesh, cl, gc, vcbc);
    Field<SymTenzor<double,D>,D,M,Cell> symtCell(mesh, cl, gc, stcbc);
    Field<Tenzor<double,D>,D,M,Cell> tensorCell(mesh, cl, gc, tcbc);
    Field<double,D,M,Cell> scalarCell(mesh, cl, gc, scbc);

    //---------------------------------------------------------------------------

    //---------------------------------------------------------------------------
    // Try out Divergence Vektor/Vert -> Scalar/Cell:
    // Assign values into the vert-centered Field<Vektor>:
    assign(vectorVert, mesh.getVertexPositionField(vectorVert));
    scalarCell = Div(vectorVert, scalarCell);
    // The value should be 3.0 for all elements; test this:
124
    EXPECT_NEAR(std::abs(sum(scalarCell)/totcells), 1.0*D, roundOffError);
125 126 127 128 129 130 131 132 133 134 135
    //---------------------------------------------------------------------------

    // --------------------------------------------------------------------------
    // Try out Gradient Scalar/Cell -> Vektor/Vert:

    // Use mesh object and vectorVert and scalarCell Field's constructed above.
    vectorCell = mesh.getCellPositionField(vectorCell);
    vectorCell -= mesh.get_origin();
    // Assign positive-sloping linear ramp values into the cell-centered
    // Field<scalar>:
    scalarCell = 0.0;
136
    for (d = 0; d < D; d++) scalarCell[cells] += vectorCell[cells](d);
137 138 139 140 141 142
    // Now take the gradient:
    vectorVert = Grad(scalarCell, vectorVert);
    // The value should be (1.0,1.0,1.0) for all elements one at least one
    // removed from the last-physical-layer elements. Last-physical-layer
    // elements will be different because the BC available in IPPL don't really
    // do the kind of linear extrapolation appropriate for the needs here:
143 144
    Vektor<double,D> unit;
    for (d = 0; d < D; d++) unit[d] = 1.0;
145 146 147 148 149
    Vektor<double,D> sumVectorVert;
    // Use temporary, smaller BareField as a reduced-by-two vector Field to hold
    // only the boundary-exclusive elements (needed because of limitations of
    // IPPL reductions ops):
    NDIndex<D> bev;
150
    for (d = 0; d < D; d++) bev[d] = Index(1,nverts[d]-2,1);
151 152 153 154
    FieldLayout<D> templayout(bev);
    BareField<Vektor<double,D>,D> temp(templayout);
    temp[bev] = vectorVert[bev];
    sumVectorVert = sum(temp);
155 156
    unsigned totred=1;
    for (d = 0; d < D; d++) totred *= nverts[d] - 2;
157 158 159 160
    sumVectorVert /= totred;
    Vektor<double,D> diffVectorVert;
    diffVectorVert = sumVectorVert - unit;
    double magDiffVectorVert = 0.0;
161
    for (d = 0; d < D; d++) magDiffVectorVert += diffVectorVert(d)*diffVectorVert(d);
162
    magDiffVectorVert = sqrt(magDiffVectorVert);
163
    EXPECT_NEAR(std::abs(magDiffVectorVert), 0, roundOffError);
164 165 166 167 168 169 170 171 172 173 174
    //---------------------------------------------------------------------------

    // --------------------------------------------------------------------------
    // Try out Gradient Scalar/Cell -> Vektor/Cell:

    // Use mesh object and vectorVert and scalarCell Field's constructed above.
    vectorCell = mesh.getCellPositionField(vectorCell);
    vectorCell -= mesh.get_origin();
    // Assign positive-sloping linear ramp values into the cell-centered
    // Field<scalar>:
    scalarCell = 0.0;
175
    for (d = 0; d < D; d++) scalarCell[cells] += vectorCell[cells](d);
176 177 178 179 180 181
    // Now take the gradient:
    vectorCell = Grad(scalarCell, vectorCell);
    // The value should be (1.0,1.0,1.0) for all elements one at least one
    // removed from the last-physical-layer elements. Last-physical-layer
    // elements will be different because the BC available in IPPL don't really
    // do the kind of linear extrapolation appropriate for the needs here:
182
    for (d = 0; d < D; d++) unit[d] = 1.0;
183 184 185 186 187
    Vektor<double,D> sumVectorCell;
    // Use temporary, smaller BareField as a reduced-by-two vector Field to hold
    // only the boundary-exclusive elements (needed because of limitations of
    // IPPL reductions ops):
    NDIndex<D> bec;
188
    for (d = 0; d < D; d++) bec[d] = Index(1,ncells[d]-2,1);
189 190 191 192
    FieldLayout<D> templayout2(bec);
    BareField<Vektor<double,D>,D> temp2(templayout);
    temp2[bec] = vectorCell[bec];
    sumVectorCell = sum(temp2);
193 194
    unsigned totredc=1;
    for (d = 0; d < D; d++) totredc *= ncells[d] - 2;
195 196 197 198
    sumVectorCell /= totredc;
    Vektor<double,D> diffVectorCell;
    diffVectorCell = sumVectorCell - unit;
    double magDiffVectorCell = 0.0;
199
    for (d = 0; d < D; d++) magDiffVectorCell += diffVectorCell(d)*diffVectorCell(d);
200
    magDiffVectorCell = sqrt(magDiffVectorCell);
201
    EXPECT_NEAR(std::abs(magDiffVectorCell), 0, roundOffError);
202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
    //---------------------------------------------------------------------------

    //---------------------------------------------------------------------------
    // Try out Divergence SymTenzor/Cell -> Vektor/Vert:

    // Use CenteredFieldLayout's from above object to construct SymTenzor Field:
    // Assign values into the cell-centered Field<SymTenzor>; use values from
    // cell-centered scalar Field scalarCell set up above:
    SymTenzor<double,D> unitSymTenzor = 1.0;
    symtCell = unitSymTenzor*scalarCell;
    // Now take the divergence:
    vectorVert = Div(symtCell, vectorVert);
    // The value should be (D,D,D,....) for all elements; test this:
    // Use temporary, smaller BareField as a reduced-by-two symtensor Field to
    // hold only the boundary-exclusive elements (needed because of limitations
    // of IPPL reductions ops):
    temp[bev] = vectorVert[bev];
    sumVectorVert = sum(temp);
    sumVectorVert /= totred;
221 222
    Vektor<double,D> deesVector;
    for (d = 0; d < D; d++) deesVector(d) = 1.0*D;
223 224
    diffVectorVert = sumVectorVert - deesVector;
    magDiffVectorVert = 0.0;
225
    for (d = 0; d < D; d++) magDiffVectorVert += diffVectorVert(d)*diffVectorVert(d);
226
    magDiffVectorVert = sqrt(magDiffVectorVert);
227
    EXPECT_NEAR(std::abs(magDiffVectorCell), 0, roundOffError);
228 229 230 231 232 233 234 235 236 237 238 239 240
    //---------------------------------------------------------------------------

    // --------------------------------------------------------------------------
    // Try out Gradient Vektor/Vert -> Tenzor/Cell:

    // Set up input values in Vektor/Vert field:
    vectorVert = mesh.getVertexPositionField(vectorVert);
    // Now take the gradient:
    tensorCell = Grad(vectorVert, tensorCell);
    // Since this is the gradient of the position vector (x*x_hat + y* y_hat +
    // z*z_hat), the result should be the identity tensor (NRL Plasma Formulary
    // Vector Identities section):
    Tenzor<double,D> identityTensor = 0.0;
241
    for (d = 0; d < D; d++) identityTensor(d,d) = 1.0;
242 243 244 245 246
    Tenzor<double,D> sumTensorCell = sum(tensorCell);
    sumTensorCell /= totcells;
    Tenzor<double,D> diffTensorCell;
    diffTensorCell = sumTensorCell - identityTensor;
    double magDiffTensorCell = 0.0;
247 248
    for (d = 0; d < D; d++) {
        for (unsigned int d2 = 0; d2 < D; d2++) {
249 250 251
            magDiffTensorCell += diffTensorCell(d,d2)*diffTensorCell(d,d2);
        }
    }
252
    EXPECT_NEAR(std::abs(magDiffTensorCell), 0, roundOffError);
253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
    //---------------------------------------------------------------------------

    //---------------------------------------------------------------------------

    /* THIS TEST DOES NOT COMPILE

    // Test Average() functions:
    // Scalar Field, Cell-Centered
    Field<double,D,Cartesian<D>,Cell> C(mesh, cl, gc);
    C = 1.0;
    // Scalar weight Field, Cell-Centered
    Field<double,D,Cartesian<D>,Cell> wC(mesh, cl, gc);
    wC = 2.0;
    // Scalar Field, Vert-Centered
    Field<double,D,Cartesian<D>,Vert> V(mesh, vl, gc);
    V = 1.0;
    // Scalar weight Field, Vert-Centered
    Field<double,D,Cartesian<D>,Vert> wV(mesh, vl, gc);
    wV = 2.0;
    // Field's to hold weighted averages:
    Field<double,D,Cartesian<D>,Cell> avgToC(mesh, cl, gc);
    Field<double,D,Cartesian<D>,Vert> avgToV(mesh, vl, gc);

    assign(avgToV, Average(C, wC, avgToV));
    assign(avgToC, Average(V, wV, avgToC));

    // Weighted average from Cell to Vert:
    // ada  does not work assign(avgToV, Average(C, wC, avgToV));
    // Weighted average from Vert to Cell:
    // ada dones not work assign(avgToC, Average(V, wV, avgToC));
    // Check results:
    if (sum(avgToV) != totverts) {
    testmsg << "avgToV values wrong" << endl;
    testmsg << "sum(avgToV) = " << sum(avgToV) << " ; totverts = " << totverts
    << endl;
    // passed = false;
    }
    if (sum(avgToC) != totcells) {
    testmsg << "avgToC values wrong" << endl;
    testmsg << "sum(avgToC) = " << sum(avgToC) << " ; totcells = " << totcells
    << endl;
    // passed = false;
    }

    */

    //---------------------------------------------------------------------------

    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    //---------------------------------------------------------------------------
    // Some accessor function tests:
    double theVolume, theVolume2, theVolume3;
    NDIndex<D> ndi;
    ndi[0] = Index(2,2,1);
    ndi[1] = Index(2,2,1);
    ndi[2] = Index(2,2,1);
    theVolume = mesh.getCellVolume(ndi);
    ndi[0] = Index(0,2,1);
    ndi[1] = Index(0,2,1);
    ndi[2] = Index(0,2,1);
    theVolume2 = mesh.getCellRangeVolume(ndi);
    EXPECT_NEAR(theVolume2, (6.0*12.0*18.0), roundOffError);

    ndi[0] = Index(0,3,1);
    ndi[1] = Index(0,3,1);
    ndi[2] = Index(0,3,1);
    theVolume3 = mesh.getVertRangeVolume(ndi);
    EXPECT_NEAR(theVolume3, theVolume2, roundOffError);
    //---------------------------------------------------------------------------
    Field<double,D,Cartesian<D>,Cell> theVolumes(mesh, cl);
    mesh.getCellVolumeField(theVolumes);
    EXPECT_NEAR((sum(theVolumes)/totcells), theVolume, roundOffError);
    //---------------------------------------------------------------------------
    Vektor<double,D> v;
    v(0) = 1.5; v(1) = 4.5; v(2) = 9.5;
    ndi = mesh.getNearestVertex(v);
    // nearest vertex should be (1,1,1)
331
    for (unsigned int i = 0; i < D; i++) {
332 333 334 335 336 337 338
        EXPECT_EQ((int)ndi[0].first(),  1);
        EXPECT_EQ((int)ndi[0].length(), 1);
    }
    //---------------------------------------------------------------------------
    Vektor<double,D> v1;
    v1 = mesh.getVertexPosition(ndi);
    v(0) = 2.0; v(1) = 4.0; v(2) = 12.0; // Correct value
339
    for (unsigned int i = 0; i < D; i++) {
340 341 342 343 344 345 346 347 348 349 350 351
        EXPECT_NEAR(v1(i), v(i), roundOffError);
    }
    //---------------------------------------------------------------------------
    CenteredFieldLayout<D,Cartesian<D>,Vert>
        clVert(mesh);
    Field<Vektor<double,D>,D,Cartesian<D>,Vert>
        thePositions(clVert);
    mesh.getVertexPositionField(thePositions);
    //---------------------------------------------------------------------------
    v = mesh.getDeltaVertex(ndi);
    Vektor<double,D> vcorrect;
    vcorrect(0) = 2.0; vcorrect(1) = 4.0; vcorrect(2) = 9.0;
352
    for (unsigned int i = 0; i < D; i++) {
353 354 355
        EXPECT_NEAR(vcorrect(i), v(i), roundOffError);
    }
}