Test initial remapping between grids and parcels

Use quadrature points for finding the connectivity between grids and parcels instead of MLPACK library.

Author: dongli

external/geomtk

@@ -19,10 +19,6 @@ class AdvectionManager
     double radialMixing;    // control the radial mixing degree
     double lateralMixing;   // control the lateral mixing degree
     double restoreFactor;   // control the base density restore degree
-
-    Tree *cellTree;
-    mat cellCoords;
-    vector<size_t> cellCoordsMap;
 public:
     AdvectionManager();
     virtual ~AdvectionManager();

src/AdvectionManager.cpp

@@ -7,8 +7,6 @@ Cartesian3DTest() {
 
 Cartesian3DTest::
 ~Cartesian3DTest() {
-    delete _mesh;
-    delete _domain;
 }
 
 void Cartesian3DTest::
@@ -33,7 +31,9 @@ init(const ConfigManager &configManager, AdvectionManager &advectionManager) {
     // Initialize velocity field.
     velocityField.create(*_mesh, true, true);
     io.file(dataIdx).addField("double", FULL_DIMENSION,
-                              {&velocityField(0), &velocityField(1)});
+                              {&velocityField(0),
+                               &velocityField(1),
+                               &velocityField(2)});
     // Initialize advection manager.
     advectionManager.init(configManager, *_mesh);
     // Initialize density fields for input and output.
@@ -54,33 +54,42 @@ init(const ConfigManager &configManager, AdvectionManager &advectionManager) {
 
 void Cartesian3DTest::
 setInitialCondition(AdvectionManager &advectionManager) {
+    TimeLevelIndex<2> timeIdx;
     io.open(dataIdx);
+    io.input<double>(dataIdx, timeIdx, {&velocityField(0),
+                                        &velocityField(1),
+                                        &velocityField(2)});
     double *q = new double[_mesh->totalNumGrid(CENTER)*densities.size()];
-    int j = 0;\
-    for (int t = 0; t < densities.size(); ++t) {
+    int j = 0;
+    for (arma::uword t = 0; t < densities.size(); ++t) {
         io.input<double>(dataIdx, {&densities[t]});
-        for (int i = 0; i < _mesh->totalNumGrid(CENTER); ++i) {
+        for (arma::uword i = 0; i < _mesh->totalNumGrid(CENTER); ++i) {
             q[j++] = densities[t](i);
         }
     }
-    TimeLevelIndex<2> timeIdx;
     advectionManager.input(timeIdx, q);
+    io.close(dataIdx);
 } // setInitialCondition
 
 void Cartesian3DTest::
-advanceDynamics(AdvectionManager &advectionManager) {
-    // Read in WRF-LES 3D flow.
+advanceDynamics(const TimeLevelIndex<2> &timeIdx,
+                AdvectionManager &advectionManager) {
+    io.open(dataIdx);
+    io.input<double>(dataIdx, timeIdx, {&velocityField(0),
+                                        &velocityField(1),
+                                        &velocityField(2)});
+    io.close(dataIdx);
+    Interface::advanceDynamics(timeIdx, advectionManager);
 } // advanceDynamics
 
 void Cartesian3DTest::
 output(const TimeLevelIndex<2> &timeIdx, AdvectionManager &advectionManager) {
     io.create(outputIdx);
-    for (int t = 0; t < densities.size(); ++t) {
-        for (int i = 0; i < _mesh->totalNumGrid(CENTER); ++i) {
-            std::cout << advectionManager.density(timeIdx, t, i) << std::endl;
+    for (arma::uword t = 0; t < densities.size(); ++t) {
+        for (arma::uword i = 0; i < _mesh->totalNumGrid(CENTER); ++i) {
             densities[t](i) = advectionManager.density(timeIdx, t, i);
         }
-        io.output<double>(outputIdx, timeIdx, {&densities[t]});
+        io.output<double>(outputIdx, {&densities[t]});
     }
     io.close(outputIdx);
 } // output

src/AdvectionManager.h

@@ -7,9 +7,10 @@
 
 class Cartesian3DTest
 : public geomtk::AdvectionTestInterface<2, Domain, Mesh, Field<double>, VelocityField, IOManager> {
-    int dataIdx;
+    arma::uword dataIdx;
 public:
     typedef geomtk::AdvectionManagerInterface<2, Domain, Mesh, VelocityField> AdvectionManager;
+    typedef geomtk::AdvectionTestInterface<2, Domain, Mesh, Field<double>, VelocityField, IOManager> Interface;
 
     Cartesian3DTest();
     virtual ~Cartesian3DTest();
@@ -21,7 +22,8 @@ class Cartesian3DTest
     setInitialCondition(AdvectionManager &advectionManager);
 
     virtual void
-    advanceDynamics(AdvectionManager &advectionManager);
+    advanceDynamics(const TimeLevelIndex<2> &timeIdx,
+                    AdvectionManager &advectionManager);
 
     virtual void
     output(const TimeLevelIndex<2> &timeIdx,

src/Cartesian3DTest.cpp

@@ -9,7 +9,7 @@ MeshAdaptor() {
 
 MeshAdaptor::
 ~MeshAdaptor() {
-    for (int t = 0; t < _masses.size(); ++t) {
+    for (uword t = 0; t < _masses.size(); ++t) {
         delete _masses[t];
     }
 }
@@ -19,6 +19,7 @@ init(const Mesh &mesh) {
     this->mesh = &mesh;
     _numConnectedParcel.resize(mesh.totalNumGrid(CENTER));
     _connectedParcels.resize(mesh.totalNumGrid(CENTER));
+    _numContainedQuadraturePoint.resize(mesh.totalNumGrid(CENTER));
     _remapWeights.resize(mesh.totalNumGrid(CENTER));
     _numContainedParcel.resize(mesh.totalNumGrid(CENTER));
     _containedParcels.resize(mesh.totalNumGrid(CENTER));
@@ -32,43 +33,64 @@ addTracer(const string &name, const string &unit, const string &comment) {
 
 void MeshAdaptor::
 resetTracers(const TimeLevelIndex<2> &timeIdx) {
-    for (int t = 0; t < _masses.size(); ++t) {
-        for (int i = 0; i < mesh->totalNumGrid(CENTER, mesh->domain().numDim()); ++i) {
-            (*_masses[t])(timeIdx, i);
+    for (uword t = 0; t < _masses.size(); ++t) {
+        for (uword i = 0; i < mesh->totalNumGrid(CENTER, mesh->domain().numDim()); ++i) {
+            (*_masses[t])(timeIdx, i) = 0;
         }
     }
 }
 
 void MeshAdaptor::
 connectParcel(int cellIdx, Parcel *parcel, double weight) {
-#ifndef NDEBUG
-    for (int i = 0; i < _numConnectedParcel[cellIdx]; ++i) {
+    for (uword i = 0; i < _numConnectedParcel[cellIdx]; ++i) {
         if (_connectedParcels[cellIdx][i] == parcel) {
-            REPORT_ERROR("Parcel (ID = " << parcel->ID() <<
-                         ") has already been connected!");
+            _numContainedQuadraturePoint[cellIdx][i]++;
+            _remapWeights[cellIdx][i] += weight;
+            return;
         }
     }
-#endif
     if (_numConnectedParcel[cellIdx] == _connectedParcels[cellIdx].size()) {
         _connectedParcels[cellIdx].push_back(parcel);
+        _numContainedQuadraturePoint[cellIdx].push_back(1);
         _remapWeights[cellIdx].push_back(weight);
     } else {
         _connectedParcels[cellIdx][_numConnectedParcel[cellIdx]] = parcel;
+        _numContainedQuadraturePoint[cellIdx][_numConnectedParcel[cellIdx]] = 1;
         _remapWeights[cellIdx][_numConnectedParcel[cellIdx]] = weight;
     }
     _numConnectedParcel[cellIdx]++;
 }
 
 void MeshAdaptor::
 resetConnectedParcels() {
-    for (int i = 0; i < mesh->totalNumGrid(CENTER, mesh->domain().numDim()); ++i) {
+    for (uword i = 0; i < mesh->totalNumGrid(CENTER, mesh->domain().numDim()); ++i) {
         _numConnectedParcel[i] = 0;
     }
 }
 
+uword MeshAdaptor::
+numContainedQuadraturePoint(int cellIdx, Parcel *parcel) const {
+    for (uword i = 0; i < _numConnectedParcel[cellIdx]; ++i) {
+        if (_connectedParcels[cellIdx][i] == parcel) {
+            return _numContainedQuadraturePoint[cellIdx][i];
+        }
+    }
+    REPORT_ERROR("Parcel is not connected!");
+}
+
 double MeshAdaptor::
 remapWeight(int cellIdx, Parcel *parcel) const {
-    for (int i = 0; i < _numConnectedParcel[cellIdx]; ++i) {
+    for (uword i = 0; i < _numConnectedParcel[cellIdx]; ++i) {
+        if (_connectedParcels[cellIdx][i] == parcel) {
+            return _remapWeights[cellIdx][i];
+        }
+    }
+    REPORT_ERROR("Parcel is not connected!");
+}
+
+double& MeshAdaptor::
+remapWeight(int cellIdx, Parcel *parcel) {
+    for (uword i = 0; i < _numConnectedParcel[cellIdx]; ++i) {
         if (_connectedParcels[cellIdx][i] == parcel) {
             return _remapWeights[cellIdx][i];
         }
@@ -79,7 +101,7 @@ remapWeight(int cellIdx, Parcel *parcel) const {
 void MeshAdaptor::
 containParcel(int cellIdx, Parcel *parcel) {
 #ifndef NDEBUG
-    for (int i = 0; i < _numContainedParcel[cellIdx]; ++i) {
+    for (uword i = 0; i < _numContainedParcel[cellIdx]; ++i) {
         if (_containedParcels[cellIdx][i] == parcel) {
             REPORT_ERROR("Parcel (ID = " << parcel->ID() <<
                          ") has already been contained!");
@@ -91,13 +113,13 @@ containParcel(int cellIdx, Parcel *parcel) {
     } else {
         _containedParcels[cellIdx][_numContainedParcel[cellIdx]] = parcel;
     }
-    parcel->hostCellIdx() = cellIdx;
+    parcel->hostCellIndex() = cellIdx;
     _numContainedParcel[cellIdx]++;
 }
 
 void MeshAdaptor::
 resetContainedParcels() {
-    for (int i = 0; i < mesh->totalNumGrid(CENTER, mesh->domain().numDim()); ++i) {
+    for (uword i = 0; i < mesh->totalNumGrid(CENTER, mesh->domain().numDim()); ++i) {
         _numContainedParcel[i] = 0;
     }
 }

src/Cartesian3DTest.h

@@ -8,10 +8,11 @@ namespace lasm {
 class MeshAdaptor {
     const Mesh *mesh;
     vector<Field<double, 2>*> _masses;
-    vector<int> _numConnectedParcel;
+    vector<uword> _numConnectedParcel;
     vector<vector<Parcel*> > _connectedParcels;
+    vector<vector<uword> > _numContainedQuadraturePoint;
     vector<vector<double> > _remapWeights;
-    vector<int> _numContainedParcel;
+    vector<uword> _numContainedParcel;
     vector<vector<Parcel*> > _containedParcels;
 public:
     MeshAdaptor();
@@ -52,10 +53,16 @@ class MeshAdaptor {
     void
     resetConnectedParcels();
 
+    uword
+    numContainedQuadraturePoint(int cellIdx, Parcel *parcel) const;
+
     double
     remapWeight(int cellIdx, Parcel *parcel) const;
 
-    int
+    double&
+    remapWeight(int cellIdx, Parcel *parcel);
+
+    uword
     numConnectedParcel(int cellIdx) const {
         return _numConnectedParcel[cellIdx];
     }
@@ -70,7 +77,7 @@ class MeshAdaptor {
     void
     resetContainedParcels();
 
-    int
+    uword
     numContainedParcel(int cellIdx) const {
         return _numContainedParcel[cellIdx];
     }

src/MeshAdaptor.cpp

@@ -36,7 +36,7 @@ init(int ID) {
         _detH.level(l) = -999.0;
         _invH.level(l).set_size(domain->numDim(), domain->numDim());
         _shapeSize.level(l).set_size(domain->numDim());
-        _meshIdx.level(l).setNumDim(domain->numDim());
+        _meshIndex.level(l).init(domain->numDim());
     }
     _U.set_size(domain->numDim(), domain->numDim());
     _V.set_size(domain->numDim(), domain->numDim());
@@ -53,7 +53,7 @@ updateDeformMatrix(const TimeLevelIndex<2> &timeIdx) {
     mat &H = _H.level(timeIdx);
     const field<BodyCoord> &y = SkeletonPoints::bodyCoords();
     const field<SpaceCoord> &x = _skeletonPoints->localSpaceCoords(timeIdx);
-    const field<MeshIndex> &meshIdxs = _skeletonPoints->meshIdxs(timeIdx);
+    const field<MeshIndex> &meshIndexs = _skeletonPoints->meshIndexs(timeIdx);
     if (domain->numDim() == 2) {
         // Calculate the elements of four matrices.
         double h11_1 = x[0](0)/y[0](0); double h21_1 = x[0](1)/y[0](0);
@@ -69,22 +69,22 @@ updateDeformMatrix(const TimeLevelIndex<2> &timeIdx) {
         double h21_1 = 0, h21_3 = 0, h22_2 = 0, h22_4 = 0, h23_5 = 0, h23_6 = 0;
         double h31_1 = 0, h31_3 = 0, h32_2 = 0, h32_4 = 0, h33_5 = 0, h33_6 = 0;
         double w1 = 0, w2 = 0, w3 = 0, w4 = 0, w5 = 0, w6 = 0;
-        if (meshIdxs[0].isValid()) {
+        if (meshIndexs[0].isValid()) {
             w1 = 1; h11_1 = x[0](0)/y[0](0); h21_1 = x[0](1)/y[0](0); h31_1 = x[0](2)/y[0](0);
         }
-        if (meshIdxs[1].isValid()) {
+        if (meshIndexs[1].isValid()) {
             w2 = 1; h12_2 = x[1](0)/y[1](1); h22_2 = x[1](1)/y[1](1); h32_2 = x[1](2)/y[1](1);
         }
-        if (meshIdxs[2].isValid()) {
+        if (meshIndexs[2].isValid()) {
             w3 = 1; h11_3 = x[2](0)/y[2](0); h21_3 = x[2](1)/y[2](0); h31_3 = x[2](2)/y[2](0);
         }
-        if (meshIdxs[3].isValid()) {
+        if (meshIndexs[3].isValid()) {
             w4 = 1; h12_4 = x[3](0)/y[3](1); h22_4 = x[3](1)/y[3](1); h32_4 = x[3](2)/y[3](1);
         }
-        if (meshIdxs[4].isValid()) {
+        if (meshIndexs[4].isValid()) {
             w5 = 1; h13_5 = x[4](0)/y[4](2); h23_5 = x[4](1)/y[4](2); h33_5 = x[4](2)/y[4](2);
         }
-        if (meshIdxs[5].isValid()) {
+        if (meshIndexs[5].isValid()) {
             w6 = 1; h13_6 = x[5](0)/y[5](2); h23_6 = x[5](1)/y[5](2); h33_6 = x[5](2)/y[5](2);
         }
         // The final matrix is the combination of the six.
@@ -118,30 +118,30 @@ updateDeformMatrix(const TimeLevelIndex<2> &timeIdx, const vec &S) {
     longAxisVertexY() = _invH.level(timeIdx)*_H.level(timeIdx)*_V.col(0);
     calcSpaceCoord(timeIdx, longAxisVertexY, longAxisVertexX);
     _filamentDegree.level(timeIdx) = _S[0]/_S.min();
-}
+} // updateDeformMatrix
 
 void Parcel::
 resetSkeletonPoints(const TimeLevelIndex<2> &timeIdx, const Mesh &mesh) {
     const field<BodyCoord> &y = SkeletonPoints::bodyCoords();
     field<SpaceCoord> &x = _skeletonPoints->spaceCoords(timeIdx);
-    field<MeshIndex> &meshIdxs = _skeletonPoints->meshIdxs(timeIdx);
-    for (int i = 0; i < x.size(); ++i) {
+    field<MeshIndex> &meshIndexs = _skeletonPoints->meshIndexs(timeIdx);
+    for (uword i = 0; i < x.size(); ++i) {
         calcSpaceCoord(timeIdx, y[i], x[i]);
-        meshIdxs[i].reset();
-        meshIdxs[i].locate(mesh, x[i]);
+        meshIndexs[i].reset();
+        meshIndexs[i].locate(mesh, x[i]);
     }
 } // resetSkeletonPoints
 
 void Parcel::
 updateVolume(const TimeLevelIndex<2> &timeIdx, double volume) {
     _detH.level(timeIdx) = volume;
-}
+} // updateVolume
 
 void Parcel::
 updateShapeSize(const TimeLevelIndex<2> &timeIdx) {
     BodyCoord y(domain->numDim());
     SpaceCoord x(domain->numDim());
-    for (int m = 0; m < domain->numDim(); ++m) {
+    for (uword m = 0; m < domain->numDim(); ++m) {
         y() = _invH.level(timeIdx)*_H.level(timeIdx)*_V.col(m);
         calcSpaceCoord(timeIdx, y, x);
         _shapeSize.level(timeIdx)[m] = domain->calcDistance(x, _x.level(timeIdx));
@@ -171,10 +171,15 @@ shapeFunction(const TimeLevelIndex<2> &timeIdx,
 
 void Parcel::
 connectCell(int cellIdx) {
-    if (_numConnectedCell == _connectedCellIdxs.size()) {
-        _connectedCellIdxs.push_back(cellIdx);
+    for (uword i = 0; i < _numConnectedCell; ++i) {
+        if (_connectedCellIndexs[i] == cellIdx) {
+            return;
+        }
+    }
+    if (_numConnectedCell == _connectedCellIndexs.size()) {
+        _connectedCellIndexs.push_back(cellIdx);
     } else {
-        _connectedCellIdxs[_numConnectedCell] = cellIdx;
+        _connectedCellIndexs[_numConnectedCell] = cellIdx;
     }
     _numConnectedCell++;
 } // connectCell