Bugfix: simple matrix operations

src/foam/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -258,8 +258,7 @@ public:
             );
 
 
-        // Destructor
-
+        //- Destructor
         virtual ~solver()
         {}
 
@@ -418,8 +417,7 @@ public:
             );
 
 
-        // Destructor
-
+        //- Destructor
         virtual ~smoother()
         {}
 
@@ -556,8 +554,7 @@ public:
             );
 
 
-        // Destructor
-
+        //- Destructor
         virtual ~preconditioner()
         {}
 
@@ -635,8 +632,7 @@ public:
         }
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~lduMatrix();
 
 
@@ -700,6 +696,11 @@ public:
                 return (lowerPtr_);
             }
 
+            bool empty() const
+            {
+                return (!diagPtr_ && !lowerPtr_ && !upperPtr_);
+            }
+
             bool diagonal() const
             {
                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);

src/foam/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C

@@ -143,6 +143,12 @@ void Foam::lduMatrix::negate()
 
 void Foam::lduMatrix::operator+=(const lduMatrix& A)
 {
+    // Escape empty matrix
+    if (A.empty())
+    {
+        return;
+    }
+    
     if (A.diagPtr_)
     {
         diag() += A.diag();
@@ -211,6 +217,12 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A)
 
 void Foam::lduMatrix::operator-=(const lduMatrix& A)
 {
+    // Escape empty matrix
+    if (A.empty())
+    {
+        return;
+    }
+    
     if (A.diagPtr_)
     {
         diag() -= A.diag();