diff --git a/mgrid_2/getmatrix.cpp b/mgrid_2/getmatrix.cpp
index bd00fb1..17deea4 100644
--- a/mgrid_2/getmatrix.cpp
+++ b/mgrid_2/getmatrix.cpp
@@ -456,30 +456,66 @@ double FEM_Matrix::ThermalConductivity(const int subdomain)
 
 double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
 {
-    double lambda = 0.0;
+    double c = 0.0;
     switch (subdomain)
     {
         // ceramic mug
         case 1:
-            lambda = 2.0 * 1e6;
+            c = 2.0 * 1e6;
             break;
 
         // water
         case 2:
-            lambda = 4.184 * 1e6;       
+            c = 4.184 * 1e6;       
             break;
 
         // air
         case 3:
-            lambda = 1.2 * 1e3;
+            c = 1.2 * 1e3;
             break;
 
         default:
-            lambda = 1.0;
+            c = 1.0;
             break;
     }
 
-    return lambda;
+    return c;
+}
+
+void FEM_Matrix::AddMass_mult(vector<double> &f)
+{
+    cout << "\n############   FEM_Matrix::AddMass_mult ";
+    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]);
+
+    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 c = VolumetricHeatCapacity(subdomain);
+
+        //cout << subdomain << endl;
+
+        CalcElem_MasseSpecific(ia.data() + 3 * i, xc.data(), c, 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;
 }
 
 void FEM_Matrix::CalculateLaplace(vector<double> &f)
@@ -932,6 +968,27 @@ void CalcElem_Masse(int const ial[3], double const xc[], double ske[3][3])
     return;
 }
 
+void CalcElem_MasseSpecific(int const ial[3], double const xc[], double const c, 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] += c * jac / 12.0;
+    ske[0][1] += c * jac / 24.0;
+    ske[0][2] += c * jac / 24.0;
+    ske[1][0] += c * jac / 24.0;
+    ske[1][1] += c * jac / 12.0;
+    ske[1][2] += c * jac / 24.0;
+    ske[2][0] += c * jac / 24.0;
+    ske[2][1] += c * jac / 24.0;
+    ske[2][2] += c * jac / 12.0;
+
+    return;
+}
+
 // #####################################################################
 
 BisectInterpolation::BisectInterpolation()
diff --git a/mgrid_2/getmatrix.h b/mgrid_2/getmatrix.h
index 6761ef3..eb28851 100644
--- a/mgrid_2/getmatrix.h
+++ b/mgrid_2/getmatrix.h
@@ -357,6 +357,16 @@ class FEM_Matrix: public CRS_Matrix
       //  c * rho [ J / (m^3 * K) ]
        double VolumetricHeatCapacity(const int subdomain);
 
+       /**
+        * Calculates the entries of f.e. mass matrix
+        * for multiple domains with different heat capacities
+        * and load/rhs vector @p f.
+        * No memory is allocated.
+        *
+        * @param[in,out] f (preallocated) rhs/load vector
+        */
+      void AddMass_mult(std::vector<double> &f);
+
         /**
         * Calculates the entries of f.e. stiffness matrix for the Laplace operator
         * and load/rhs vector @p f.
@@ -698,6 +708,16 @@ void CalcElemSpecific(int const ial[3], double const xc[], double const lambda,
  */
 void CalcElem_Masse(int const ial[3], double const xc[], double ske[3][3]);
 
+/**
+ * Calculates the element mass matrix @p ske.
+ * of one triangular element with linear shape functions
+ * for specific volumetric heat capacity 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[out]  ske element stiffness matrix
+ */
+void CalcElem_MasseSpecific(int const ial[3], double const xc[], double const c, double ske[3][3]);
+
 /**
  * Calculates element load vector @p fe of one triangular element with linear shape functions.
  * @param[in]	ial	 node indices of the three element vertices
diff --git a/mgrid_2/main.cpp b/mgrid_2/main.cpp
index 576b0a9..973245a 100644
--- a/mgrid_2/main.cpp
+++ b/mgrid_2/main.cpp
@@ -45,17 +45,19 @@ int main(int argc, char **argv )
     // Initialize RHS
     vector<double> fv(SK.Nrows(), 0.0);    // r.h.s.
 
-    // Calculate Matrix entries
+    // Calculate stiffness matrix entries
     SK.CalculateLaplaceMult(fv);                       // matrix
     //SK.Debug();
 
+    // Add mass matrix entries
+    SK.AddMass_mult(fv);
+
     // 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.CalculateRHS(fv, [](double x, double y) {return 0;});
     //SK.CheckRowSum();
     SK.CheckMatrix();
 
+
     // Initialize temperature
     vector<double> uv(SK.Nrows(), 0.0);    // temperature
     mesh_c.Init_Solution_mult(uv, 0, [](double x, double y) -> double { return 18; });    // mug