From 0bb9c726053232f8abdb9680f5be8ca1685eb693 Mon Sep 17 00:00:00 2001
From: orlandini <francisco.orlandini@gmail.com>
Date: Tue, 28 May 2024 17:45:53 -0300
Subject: [PATCH] feat(elmat): adds option for user allocated mat

---
 Mesh/TPZElementMatrixT.cpp | 18 +++++++++++-------
 Mesh/TPZElementMatrixT.h   |  2 ++
 2 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/Mesh/TPZElementMatrixT.cpp b/Mesh/TPZElementMatrixT.cpp
index cd102d82e..7909ca4a3 100644
--- a/Mesh/TPZElementMatrixT.cpp
+++ b/Mesh/TPZElementMatrixT.cpp
@@ -249,6 +249,10 @@ void TPZElementMatrixT<TVar>::ApplyConstraints(){
 	
 	int numnodes_processed = 0;
 	int current_order = 0;
+  TPZFNMatrix<150,TVar> localmat;
+
+  TPZFMatrix<TVar> *fulldepmat = fUserAllocMat ? fUserAllocMat : &localmat;
+  
 	while(numnodes_processed < totalnodes) {
 		int in;
 		for(in=0; in<totalnodes; in++) {
@@ -268,7 +272,7 @@ void TPZElementMatrixT<TVar>::ApplyConstraints(){
 			// loop over the nodes from which dfn depends
 			TPZConnect::TPZDependBase *dep = dfn->FirstDepend();
 
-      TPZFNMatrix<150,TVar> fulldepmat;
+      
 			while(dep) {
         auto dept = dynamic_cast<TPZConnect::TPZDepend<TVar>*>(dep);
         if(!dept){DebugStop();}
@@ -304,7 +308,7 @@ void TPZElementMatrixT<TVar>::ApplyConstraints(){
           const int orig_col=depmat.Cols();
           const int full_row=numstate*orig_row;
           const int full_col=numstate*orig_col;
-          fulldepmat.Redim(full_row,full_col);
+          fulldepmat->Redim(full_row,full_col);
           
           for(int ic = 0; ic < orig_col; ic++){
             const auto first_c = ic*numstate;
@@ -312,7 +316,7 @@ void TPZElementMatrixT<TVar>::ApplyConstraints(){
               const auto first_r = ir*numstate;
               const auto val = depmat.GetVal(ir,ic);
               for(int istate = 0; istate < numstate; istate++){
-                fulldepmat.PutVal(first_r+istate,first_c+istate,val);
+                fulldepmat->PutVal(first_r+istate,first_c+istate,val);
               }
             }
           }
@@ -327,12 +331,12 @@ void TPZElementMatrixT<TVar>::ApplyConstraints(){
           since there are no trial functions
         */
 
-        const auto deprows = fulldepmat.Rows();
-        const auto depcols = fulldepmat.Cols();
+        const auto deprows = fulldepmat->Rows();
+        const auto depcols = fulldepmat->Cols();
         //the window is the full matrix
 
         {
-          TPZMatrixWindow<TVar> dep_window(fulldepmat,0,0,deprows,depcols);
+          TPZMatrixWindow<TVar> dep_window(*fulldepmat,0,0,deprows,depcols);
           TPZMatrixWindow<TVar> send_window(this->fConstrMat,inpos,0,insize,this->fConstrMat.Cols());
           TPZMatrixWindow<TVar> receive_window(this->fConstrMat,deppos,0,depsize,this->fConstrMat.Cols());
           const TVar alpha{1};
@@ -346,7 +350,7 @@ void TPZElementMatrixT<TVar>::ApplyConstraints(){
         }
         if (this->fType == TPZElementMatrix::EK){
           //now we multiply it on the right side too
-          TPZMatrixWindow<TVar> dep_window(fulldepmat,0,0,deprows,depcols);
+          TPZMatrixWindow<TVar> dep_window(*fulldepmat,0,0,deprows,depcols);
           TPZMatrixWindow<TVar> send_window(this->fConstrMat,0,inpos,this->fConstrMat.Rows(),insize);
           TPZMatrixWindow<TVar> receive_window(this->fConstrMat,0,deppos,this->fConstrMat.Rows(),depsize);
           const TVar alpha{1};
diff --git a/Mesh/TPZElementMatrixT.h b/Mesh/TPZElementMatrixT.h
index 2674e85e2..1e619d3ba 100644
--- a/Mesh/TPZElementMatrixT.h
+++ b/Mesh/TPZElementMatrixT.h
@@ -75,6 +75,7 @@ struct TPZElementMatrixT : public TPZElementMatrix {
     TPZBlock & ConstrBlock() override{
         return fConstrBlock;
     }
+  void SetUserAllocMat(TPZFMatrix<TVar> *mat){fUserAllocMat=mat;}
     /** @brief Pointer to a blocked matrix object*/
 	TPZFNMatrix<1000, TVar> fMat;
     /** @brief Block structure associated with fMat*/
@@ -83,6 +84,7 @@ struct TPZElementMatrixT : public TPZElementMatrix {
 	TPZFNMatrix<1000, TVar> fConstrMat;
     /** @brief Block structure associated with fConstrMat*/
 	TPZBlock fConstrBlock;
+  TPZFMatrix<TVar> *fUserAllocMat{nullptr};
 };
 
 extern template class TPZElementMatrixT<STATE>;