added subdomain support in Mesh class, CalculateLaplaceMult implementation

mgrid_2/geom.cpp

@@ -236,6 +236,7 @@ void Mesh::Export_scicomp(std::string const &basename) const
     return;
 }
 
+
 // subject to permutation:
 //     re-sort:  _xc        
 //               _xc[2*k_new], _xc[2*k_new+1]  with k_new = po2n[k] via old(_xc);
@@ -283,9 +284,10 @@ void Mesh::Visualize(vector<double> const &v) const
 void Mesh::Visualize_matlab(vector<double> const &v) const
 {
     // define external command
-    const string exec_m("matlab -nosplash < visualize_results.m");                 // Matlab
+    //const string exec_m("matlab -nosplash < visualize_results.m");                 // Matlab
     //const string exec_m("octave --no-window-system --no-gui visualize_results.m"); // Octave
     //const string exec_m("flatpak run org.octave.Octave visualize_results.m");      // Octave (flatpak): desktop GH
+    const string exec_m("octave visualize_results.m");
 
     const string fname("uv.txt");
     Write_ascii_matlab(fname, v);
@@ -956,6 +958,45 @@ Mesh::Mesh(std::string const &fname)
     //cout << " P E R M U T E D !" << endl;
 }
 
+
+    vector<int> ElementSubdomains;
+
+
+    Mesh::Mesh(std::string const &filename, std::string const &subdomain_filename) : Mesh(filename)
+    {
+        ElementSubdomains = ReadElementSubdomains(subdomain_filename);
+    }
+
+const vector<int> Mesh::ReadElementSubdomains(string const &dname) const
+{
+    ifstream ifs(dname);
+    if (!(ifs.is_open() && ifs.good())) {
+        cerr << "ParMesh::ReadElementSubdomain: Error cannot open file " << dname << endl;
+        assert(ifs.is_open());
+    }
+
+    int const OFFSET{1};             // Matlab to C indexing
+    cout << "ASCI file  " << dname << "  opened" << endl;
+
+    // Read some mesh constants
+    int nelem;
+    ifs >> nelem;
+    cout << nelem << "  " << Nelems() << endl;
+    assert( Nelems() == nelem);
+
+    // Allocate memory
+    vector<int> t2d(nelem, -1);
+
+    // Read element mapping
+    for (int k = 0; k < nelem; ++k) {
+        int tmp;
+        ifs >> tmp;
+        t2d[k] = tmp - OFFSET;
+    }
+
+    return t2d;
+}
+
 void Mesh::ReadVertexBasedMesh(std::string const &fname)
 {
     ifstream ifs(fname);

mgrid_2/geom.h

@@ -45,6 +45,16 @@ public:
     */
     explicit Mesh(std::string const &fname);
 
+    /**
+     * Reads mesh data plus subdomain data from a binary file.
+     *
+     * File format, see ascii_write_mesh.m
+     *
+     * @param[in] filename file name
+     * @param[in] subdomain_filename subdomain file name
+    */
+    explicit Mesh(std::string const &filename, std::string const &subdomain_filename);
+
     /**
      * Reads mesh data from a binary file.
      *
@@ -63,6 +73,8 @@ public:
         return _nelem;
     }
 
+    
+
     /**
      * Global number of vertices for each finite element.
      * @return number of vertices per element.
@@ -422,6 +434,7 @@ public:
      */
     void liftToQuadratic();    
 
+
 protected:
     //public:
     void SetNelem(int nelem)
@@ -522,7 +535,19 @@ public:
      */    
     [[nodiscard]] bool checkObtuseAngles() const;
 
-private:
+
+    std::vector<int> ElementSubdomains;
+
+    /**
+     * Reads the global triangle to subdomain mapping.
+     *
+     * @param[in] dname file name
+     *
+     * @see ascii_write_subdomains.m for the file format
+    */
+    [[nodiscard]] const std::vector<int> ReadElementSubdomains(std::string const &dname) const;
+
+
     /**
      * Calculates the largest inner angle in element @p idx.
      * 

mgrid_2/getmatrix.cpp

@@ -383,6 +383,87 @@ void FEM_Matrix::Derive_Matrix_Pattern_slow()
     return;
 }
 
+void FEM_Matrix::CalculateLaplaceMult(vector<double> &f)
+{
+    cout << "\n############   FEM_Matrix::CalculateLaplaceMult ";
+    double tstart = omp_get_wtime();                  // OpenMP
+    assert(_mesh.NdofsElement() == 3);               // only for triangular, linear elements
+    //cout << _nnz << " vs. " << _id[_nrows] << "  " << _nrows<< endl;
+    assert(_nnz == _id[_nrows]);
+
+    for (int k = 0; k < _nrows; ++k) {
+        _sk[k] = 0.0;
+    }
+    for (int k = 0; k < _nrows; ++k) {
+        f[k] = 0.0;
+    }
+
+    double ske[3][3], fe[3];
+    //  Loop over all elements
+    auto const nelem = _mesh.Nelems();
+    auto const &ia   = _mesh.GetConnectivity();
+    auto const &xc   = _mesh.GetCoords();
+
+    const vector<int> sd_vec = _mesh.ElementSubdomains;
+
+    #pragma omp parallel for private(ske,fe)
+    for (int i = 0; i < nelem; ++i) {
+
+        auto subdomain = sd_vec[i];
+        double lambda = Thermal_coefficient(subdomain);
+
+        cout << subdomain << endl;
+
+        CalcElemSpecific(ia.data() + 3 * i, xc.data(), lambda, ske);
+        //AddElem(ia.data()+3 * i, ske, fe, _id.data(), _ik.data(), _sk.data(), f.data()); // GH: deprecated
+        AddElem_3(ia.data() + 3 * i, ske, fe, f);
+    }
+
+    double duration = omp_get_wtime() - tstart;             // OpenMP
+    cout << "finished in  " <<  duration  << " sec.    ########\n";  // ToDo: change to  systemclock
+    //Debug();
+
+
+
+
+    return;
+}
+
+double FEM_Matrix::Thermal_coefficient(const int subdomain)
+{
+    int matlab_sd_index = subdomain - 1;
+    double lambda = 0.0;
+
+    switch (matlab_sd_index)
+    {
+        // outside
+        case 0:
+            lambda = 1.0;
+            break;
+
+        // ceramic
+        case 1:
+            lambda = 1.0;
+            break;
+
+        // water
+        case 2:
+            lambda = 1.0;
+            break;
+
+        // air
+        case 3:
+            lambda = 1.0;
+            break;
+
+        default:
+            lambda = 1.0;
+            break;
+    }
+
+    return lambda;
+}
+
 
 void FEM_Matrix::CalculateLaplace(vector<double> &f)
 {
@@ -686,6 +767,26 @@ void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3
 }
 
 
+void CalcElemSpecific(int const ial[3], double const xc[], double const lambda, double ske[3][3])
+{
+    const int  i1  = 2 * ial[0],   i2 = 2 * ial[1],   i3 = 2 * ial[2];
+    const double x13 = xc[i3 + 0] - xc[i1 + 0],  y13 = xc[i3 + 1] - xc[i1 + 1],
+                 x21 = xc[i1 + 0] - xc[i2 + 0],  y21 = xc[i1 + 1] - xc[i2 + 1],
+                 x32 = xc[i2 + 0] - xc[i3 + 0],  y32 = xc[i2 + 1] - xc[i3 + 1];
+    const double jac = fabs(x21 * y13 - x13 * y21);
+
+    ske[0][0] = lambda * 0.5 / jac * (y32 * y32 + x32 * x32);
+    ske[0][1] = lambda * 0.5 / jac * (y13 * y32 + x13 * x32);
+    ske[0][2] = lambda * 0.5 / jac * (y21 * y32 + x21 * x32);
+    ske[1][0] = ske[0][1];
+    ske[1][1] = lambda * 0.5 / jac * (y13 * y13 + x13 * x13);
+    ske[1][2] = lambda * 0.5 / jac * (y21 * y13 + x21 * x13);
+    ske[2][0] = ske[0][2];
+    ske[2][1] = ske[1][2];
+    ske[2][2] = lambda * 0.5 / jac * (y21 * y21 + x21 * x21);
+}
+
+
 void CalcElem_RHS(int const ial[3], double const xc[], double fe[3],
              const std::function<double(double,double)> &func)
 {

mgrid_2/getmatrix.h

@@ -340,6 +340,19 @@ class FEM_Matrix: public CRS_Matrix
        void Derive_Matrix_Pattern_slow();
 
 
+       /**
+        * Calculates the entries of f.e. stiffness matrix for the Laplace operator
+        * for multiple domains with different conductivities
+        * and load/rhs vector @p f.
+        * No memory is allocated.
+        *
+        * @param[in,out] f (preallocated) rhs/load vector
+        */
+       void CalculateLaplaceMult(std::vector<double> &f);
+
+       double Thermal_coefficient(const int subdomain);
+
+
         /**
         * Calculates the entries of f.e. stiffness matrix for the Laplace operator
         * and load/rhs vector @p f.
@@ -647,6 +660,17 @@ class BisectIntDirichlet: public BisectInterpolation
  */
 void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3]);
 
+/**
+ * Calculates the element stiffness matrix @p ske 
+ * of one triangular element with linear shape functions
+ * for specific thermal conductivity in subdomain
+ * @param[in]	ial	node indices of the three element vertices
+ * @param[in]	xc	vector of node coordinates with x(2*k,2*k+1) as coordinates of node k
+ * @param[in]  lambda   thermal conductivity of element
+ * @param[out] ske	element stiffness matrix
+ */
+void CalcElemSpecific(int const ial[3], double const xc[], double const lambda, double ske[3][3]);
+
 /**
  * Calculates the element mass matrix @p ske.
  * of one triangular element with linear shape functions.

mgrid_2/main.cpp

@@ -23,9 +23,10 @@ int main(int argc, char **argv )
     int nrefine = 0;
     if (argc > 1)  nrefine = atoi(argv[1]);
 
+
+    // generating the mesh
+    Mesh const mesh_c("../generate_mesh/coffee_cup.txt", "../generate_mesh/coffee_cup_sd.txt");
     //Mesh const mesh_c("square_tiny.txt");
-    Mesh const mesh_c("square_100.txt");
-    //Mesh const mesh_c("square.txt");
     bool ba = mesh_c.checkObtuseAngles();
     if (ba) cout << "mesh corrected" << endl;
 
@@ -34,33 +35,36 @@ int main(int argc, char **argv )
     //mesh.Debug();
     //mesh.DebugEdgeBased();
 
+
+    // Initializing FEM matrix !pattern! (only zero entries now)
     FEM_Matrix SK(mesh);                   // CRS matrix
     //SK.writeBinary("sparseMatrix.bin");
     //SK.Debug();
 
-    vector<double> uv(SK.Nrows(), 0.0);    // temperature
+
+    // Initialize RHS
     vector<double> fv(SK.Nrows(), 0.0);    // r.h.s.
 
-    SK.CalculateLaplace(fv);                       // matrix
-    SK.CalculateRHS(fv, [](double x, double y) {   // rhs
-        return std::sin(M_PI * 2.5 * y) * (M_PI * M_PI * 2.5 * 2.5 * x * x - 2);
-    }
-                   );
-    //SK.CheckRowSum();
-    SK.CheckMatrix();
+    // Calculate Matrix entries
+    SK.CalculateLaplaceMult(fv);                       // matrix
     //SK.Debug();
 
-    // Two ways to initialize the vector
-    //mesh.SetValues(uv,f_zero);             // user function
-    //mesh.SetValues(uv, [](double x, double y) -> double {return 0.0*x*y;} );  // lambda function
-    //mesh.SetValues(uv, [](double x, double y) -> double {return 5e-3*(x+1)*(y+1);} );  // lambda function
-    //
-    mesh.SetValues(uv, [](double x, double y) -> double {
-        return x *x * std::sin(2.5 * M_PI * y);
-    } );
+    // Calculate RHS
+     SK.CalculateRHS(fv, [](double x, double y) {   // rhs
+         return std::sin(M_PI * 2.5 * y) * (M_PI * M_PI * 2.5 * 2.5 * x * x - 2); });
+    //SK.CheckRowSum();
+    SK.CheckMatrix();
 
+    // Initialize temperature
+    vector<double> uv(SK.Nrows(), 0.0);    // temperature
+    mesh.SetValues(uv, [](double x, double y) -> double { return 18; } ); // initial temperature of every domain 
+
+
+    // Apply BC
     SK.ApplyDirichletBC(uv, fv);
 
+
+    // Solve
     auto exact_sol(uv);
     //SK.Mult(fv,uv);
 
@@ -73,6 +77,8 @@ int main(int argc, char **argv )
     double t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
     cout << "JacobiSolve: timing in sec. : " << t_diff << endl;
 
+
+    // Calculate error and visualize
     auto [val, idx] = findLargestAbsError(exact_sol, uv, 1e+6, 100);
 
     //mesh.Visualize(getAbsError(exact_sol, uv));