Merge a7279a8 into 2a92cad

include/mesh.hpp

@@ -195,6 +195,8 @@ class mesh_t {
 
   // plotting info for generating field vtu
   int    plotNverts=0;    // number of vertices for each plot element
+  int    plotN=0;         // degree of plot interpolation
+  int    plotNq=0;        // plotNq = plotN+1
   int    plotNp=0;        // number of plot nodes per element
   int    plotNelements=0; // number of "plot elements" per element
   int    *plotEToV;       // triangulation of plot nodes
@@ -339,6 +341,8 @@ class mesh_t {
 
   virtual dfloat MinCharacteristicLength() = 0;
 
+  virtual void PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch=nullptr)=0;
+
   void RecursiveSpectralBisectionPartition();
 
   void MassMatrixApply(occa::memory& o_q, occa::memory& o_Mq);

include/mesh/mesh2D.hpp

@@ -66,6 +66,8 @@ class meshTri2D: public mesh2D {
   void CubatureNodes();
 
   void MassMatrixKernelSetup(int Nfields);
+
+  void PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch=nullptr);
 };
 
 class meshQuad2D: public mesh2D {
@@ -84,6 +86,8 @@ class meshQuad2D: public mesh2D {
   void CubatureNodes();
 
   void MassMatrixKernelSetup(int Nfields);
+
+  void PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch=nullptr);
 };
 
 #endif

include/mesh/mesh3D.hpp

@@ -69,6 +69,8 @@ class meshTri3D: public mesh3D {
   void CubatureNodes();
 
   void MassMatrixKernelSetup(int Nfields);
+
+  void PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch=nullptr);
 };
 
 class meshQuad3D: public mesh3D {
@@ -87,6 +89,8 @@ class meshQuad3D: public mesh3D {
   void CubatureNodes();
 
   void MassMatrixKernelSetup(int Nfields);
+
+  void PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch=nullptr);
 };
 
 class meshTet3D: public mesh3D {
@@ -105,6 +109,8 @@ class meshTet3D: public mesh3D {
   void CubatureNodes();
 
   void MassMatrixKernelSetup(int Nfields);
+
+  void PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch=nullptr);
 };
 
 class meshHex3D: public mesh3D {
@@ -123,6 +129,8 @@ class meshHex3D: public mesh3D {
   void CubatureNodes();
 
   void MassMatrixKernelSetup(int Nfields);
+
+  void PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch=nullptr);
 };
 
 #endif

libs/mesh/meshPlotInterpHex3D.cpp

@@ -0,0 +1,101 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "mesh.hpp"
+#include "mesh/mesh3D.hpp"
+
+//interpolate field to plotting nodes
+void meshHex3D::PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch){
+
+  dfloat *IQ, *IIQ;
+
+  bool alloc_scratch=false;
+  if (scratch==nullptr) {
+    //if not provided with a scratch space, alloc our own
+    alloc_scratch=true;
+    IQ  = (dfloat *) malloc(plotNq*Nq*Nq*sizeof(dfloat));
+    IIQ = (dfloat *) malloc(plotNq*plotNq*Nq*sizeof(dfloat));
+  } else {
+    IQ  = scratch;
+    IIQ = scratch + plotNq*Nq*Nq;
+  }
+
+  //interpolate in r
+  for(int k=0;k<Nq;++k){
+    for(int j=0;j<Nq;++j){
+      for(int i=0;i<plotNq;++i){
+        dfloat qn = 0;
+
+        for(int m=0;m<Nq;++m){
+          const int qid = m + j*Nq + k*Nq*Nq;
+          qn += plotInterp[i*Nq+m]*q[qid];
+        }
+
+        const int id = i + j*plotNq + k*Nq*plotNq;
+        IQ[id] = qn;
+      }
+    }
+  }
+
+  //interpolate in s
+  for(int k=0;k<Nq;++k){
+    for(int j=0;j<plotNq;++j){
+      for(int i=0;i<plotNq;++i){
+        dfloat qn = 0;
+
+        for(int m=0;m<Nq;++m){
+          const int qid = i + m*plotNq + k*Nq*plotNq;
+          qn += plotInterp[j*Nq+m]*IQ[qid];
+        }
+
+        const int id = i + j*plotNq + k*plotNq*plotNq;
+        IIQ[id] = qn;
+      }
+    }
+  }
+
+  //interpolate in k and write
+  for(int k=0;k<plotNq;++k){
+    for(int j=0;j<plotNq;++j){
+      for(int i=0;i<plotNq;++i){
+        dfloat qn = 0;
+
+        for(int m=0;m<Nq;++m){
+          const int qid = i + j*plotNq + m*plotNq*plotNq;
+          qn += plotInterp[k*Nq+m]*IIQ[qid];
+        }
+
+        const int id = i + j*plotNq + k*plotNq*plotNq;
+        Iq[id] = qn;
+      }
+    }
+  }
+
+  //clean up
+  if (alloc_scratch) {
+    free(IQ); free(IIQ);
+  }
+}

libs/mesh/meshPlotInterpQuad2D.cpp

@@ -0,0 +1,78 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "mesh.hpp"
+#include "mesh/mesh2D.hpp"
+
+//interpolate field to plotting nodes
+void meshQuad2D::PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch){
+
+  dfloat *IQ;
+
+  bool alloc_scratch=false;
+  if (scratch==nullptr) {
+    //if not provided with a scratch space, alloc our own
+    alloc_scratch=true;
+    IQ  = (dfloat *) malloc(plotNq*Nq*sizeof(dfloat));
+  } else {
+    IQ  = scratch;
+  }
+
+  //interpolate in r
+  for(int j=0;j<Nq;++j){
+    for(int i=0;i<plotNq;++i){
+      dfloat qn = 0;
+
+      for(int m=0;m<Nq;++m){
+        const int qid = m + j*Nq;
+        qn += plotInterp[i*Nq+m]*q[qid];
+      }
+
+      const int id = i + j*plotNq;
+      IQ[id] = qn;
+    }
+  }
+
+  //interpolate in s
+  for(int j=0;j<plotNq;++j){
+    for(int i=0;i<plotNq;++i){
+      dfloat qn = 0;
+
+      for(int m=0;m<Nq;++m){
+        const int qid = i + m*plotNq;
+        qn += plotInterp[j*Nq+m]*IQ[qid];
+      }
+
+      const int id = i + j*plotNq;
+      Iq[id] = qn;
+    }
+  }
+
+  //clean up
+  if (alloc_scratch) {
+    free(IQ);
+  }
+}

libs/mesh/meshPlotInterpQuad3D.cpp

@@ -0,0 +1,78 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "mesh.hpp"
+#include "mesh/mesh3D.hpp"
+
+//interpolate field to plotting nodes
+void meshQuad3D::PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch){
+
+  dfloat *IQ;
+
+  bool alloc_scratch=false;
+  if (scratch==nullptr) {
+    //if not provided with a scratch space, alloc our own
+    alloc_scratch=true;
+    IQ  = (dfloat *) malloc(plotNq*Nq*sizeof(dfloat));
+  } else {
+    IQ  = scratch;
+  }
+
+  //interpolate in r
+  for(int j=0;j<Nq;++j){
+    for(int i=0;i<plotNq;++i){
+      dfloat qn = 0;
+
+      for(int m=0;m<Nq;++m){
+        const int qid = m + j*Nq;
+        qn += plotInterp[i*Nq+m]*q[qid];
+      }
+
+      const int id = i + j*plotNq;
+      IQ[id] = qn;
+    }
+  }
+
+  //interpolate in s
+  for(int j=0;j<plotNq;++j){
+    for(int i=0;i<plotNq;++i){
+      dfloat qn = 0;
+
+      for(int m=0;m<Nq;++m){
+        const int qid = i + m*plotNq;
+        qn += plotInterp[j*Nq+m]*IQ[qid];
+      }
+
+      const int id = i + j*plotNq;
+      Iq[id] = qn;
+    }
+  }
+
+  //clean up
+  if (alloc_scratch) {
+    free(IQ);
+  }
+}

libs/mesh/meshPlotInterpTet3D.cpp

@@ -0,0 +1,42 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include "mesh.hpp"
+#include "mesh/mesh3D.hpp"
+
+//interpolate field to plotting nodes
+void meshTet3D::PlotInterp(const dfloat* q, dfloat* Iq, dfloat* scratch){
+
+  //interpolate
+  for(int n=0;n<plotNp;++n){
+    dfloat qn = 0;
+
+    for(int m=0;m<Np;++m){
+      qn += plotInterp[n*Np+m]*q[m];
+    }
+    Iq[n] = qn;
+  }
+}