From 00eb072ce6b19a676c55fee1432b585ed86a4b98 Mon Sep 17 00:00:00 2001
From: pvanberg <pierreolivier.vanberg@gmail.com>
Date: Wed, 6 Mar 2019 12:26:46 +0100
Subject: [PATCH] Some fixing, still not working.

---
 include/Element.h |  2 +-
 src/Element.cpp   | 40 +++++++++++++++++++++++++++++++++------
 src/Mesh.cpp      | 48 ++++++++++++++++++++++++++++++++---------------
 src/solver.cpp    |  4 ++--
 4 files changed, 70 insertions(+), 24 deletions(-)

diff --git a/include/Element.h b/include/Element.h
index e8e303a3..c0f6957a 100644
--- a/include/Element.h
+++ b/include/Element.h
@@ -57,7 +57,7 @@ class Element {
 
         bool hasNode(const int tag);
 
-        void getFlux(Eigen::MatrixXd &Flux, const Eigen::Vector3d &a);
+        void getFlux(Eigen::VectorXd &Flux, const Eigen::Vector3d &a, std::vector<Element> &elements);
 
         // Add face to element
         Element &addFace(Face face);
diff --git a/src/Element.cpp b/src/Element.cpp
index 4a316dfd..33d7698f 100644
--- a/src/Element.cpp
+++ b/src/Element.cpp
@@ -194,18 +194,46 @@ void Element::getStiffMatrix(Eigen::MatrixXd &stiffMatrix, const Eigen::Vector3d
     }
 }
 
-void Element::getFlux(Eigen::MatrixXd &Flux, const Eigen::Vector3d &a){
+void Element::getFlux(Eigen::VectorXd &Flux, const Eigen::Vector3d &a, std::vector<Element> &elements){
     // Independent of numerical flux
-    Flux = Eigen::MatrixXd::Zero(this->numBasisFcts, this->numBasisFcts);
+    Eigen::VectorXd numDataIn;
+    Eigen::VectorXd numDataOut;
+    Flux = Eigen::VectorXd::Zero(this->numBasisFcts);
     for(unsigned int i=0; i<this->numBasisFcts; ++i) {
         for (unsigned int j = 0; j < this->numBasisFcts; ++j) {
             for(Face &face: this->faces){
-                if(face.hasNode(this->nodeTags[i]) && face.hasNode(this->nodeTags[j])){
-                    Flux(i, j) += face.getFluxInt(a);
+                if(face.boundary){
+                    // Not flux at boundaries
+                    Flux(i) += 0.0;
                 }
-                else{
-                    Flux(i, j) += 0.0;
+                else {
+                    // If face is not a boundary
+                    // We get is element other side of face.
+                    int elOutTag = face.getSecondElement(this->tag);
+                    int idElOut;
+                    for(int i=0; i< elements.size(); ++i){
+                        if(elements[i].getTag() == tag)
+                            idElOut=  i;
+                    }
+                    Element elOut = elements[idElOut];
+
+                    // If Node 'i' is in Surface
+                    if(face.hasNode(this->nodeTags[i])){
+                        // Get indice of out solution
+                        for(int k=0; k<elOut.getNumNodes(); ++k) {
+                            if (this->getNodeTags()[i] == elOut.getNodeTags()[k]) {
+                                this->getData(numDataIn);
+                                elOut.getData(numDataOut);
+                                Flux(i) += ((numDataIn(i) + numDataOut(k)) / 2.)*face.getFluxInt(a);
+                            }
+                        }
+                    }
+                    else{
+                        // Shape function is zero on surface
+                        Flux(i) += 0.0;
+                    }
                 }
+
             }
         }
     }
diff --git a/src/Mesh.cpp b/src/Mesh.cpp
index ac3db5ee..1f56d0fc 100644
--- a/src/Mesh.cpp
+++ b/src/Mesh.cpp
@@ -147,9 +147,19 @@ Mesh::Mesh(std::string name) :  name(name){
             }
             if(hasFace) {
                 face.addElement(element->getTag());
+            }
+        }
+        for(auto element = std::begin(this->elements); element!=std::end(this->elements); ++element) {
+            bool hasFace = true;
+            for (auto &node : face.getNodeTags()) {
+                if(!element->hasNode(node))
+                    hasFace = false;
+            }
+            if(hasFace) {
                 element->addFace(face);
             }
         }
+
     }
     Log("Faces loaded");
 
@@ -209,21 +219,21 @@ void Mesh::getStiffMatrix(Eigen::SparseMatrix<double> &K, const Eigen::Vector3d
 
 void Mesh::getFlux(Eigen::VectorXd &F, const Eigen::Vector3d &a) {
     int offset = 0;
-    F.resize(this->elements.size()*this->elements[0].getNumNodes());
+    F.resize(this->elements.size() * this->elements[0].getNumNodes());
     // Loop over elements
-    for(Element& el : this->elements) {
+    for (Element &el : this->elements) {
+
+        Eigen::VectorXd numFlux;
+        el.getFlux(numFlux, a, this->elements);
 
-        Eigen::Vector3d numFlux;
-        Eigen::MatrixXd elF;
-        el.getFlux(elF, a);
 
         // Depends on numerical flux : (In+Out)/2
-        Eigen::VectorXd numDataIn;
+        /*Eigen::VectorXd numDataIn;
         Eigen::VectorXd numDataOut;
         for(int i=0; i<el.getNumNodes(); ++i){
             for(Face &face: el.faces) {
                 if(face.boundary){
-                    numFlux(i) = 0.0;
+                    numFlux(i) += 0.0;
                 }
                 else {
                     int elOutTag = face.getSecondElement(el.getTag());
@@ -232,20 +242,21 @@ void Mesh::getFlux(Eigen::VectorXd &F, const Eigen::Vector3d &a) {
                         if(el.getNodeTags()[i] == el.getNodeTags()[j]){
                             el.getData(numDataIn);
                             elOut.getData(numDataOut);
-                            numFlux(i) = (numDataIn(i)+numDataOut(j))/2;
+                            numFlux(i) = (numDataIn(i)+numDataOut(j))/2.;
                         }
                     }
                 }
             }
-        }
+        }*/
 
-        Eigen::VectorXd FF = elF*numFlux;
-        F(offset+0)= FF(0);
-        F(offset+1)= FF(2);
-        F(offset+1)= FF(2);
-        F(offset+1)= FF(2);
+        //Eigen::VectorXd FF = elF*numFlux;
+        F(offset + 0) = numFlux(0);
+        F(offset + 1) = numFlux(1);
+        F(offset + 2) = numFlux(2);
         offset += el.getNumNodes();
     }
+    Eigen::VectorXd data;
+    this->getData(data);
 }
 
 Element &Mesh::getElement(const int tag){
@@ -274,7 +285,14 @@ void Mesh::getData(Eigen::VectorXd &data){
     Eigen::VectorXd elData;
     for(Element& el : this->elements) {
         el.getData(elData);
-        data << data, elData;
+        if(data.size() == 0){
+            data = elData;
+        }
+        else{
+            Eigen::VectorXd temp(data.size()+elData.size());
+            temp << data, elData;
+            data = temp;
+        }
     }
 }
 
diff --git a/src/solver.cpp b/src/solver.cpp
index 18c2f0b7..deda98ec 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -53,7 +53,7 @@ namespace solver {
         Eigen::SparseMatrix<double> M_inv = solver.solve(I);
 
         // Compute stiffness (or convection) matrix
-        Eigen::Vector3d a = {1, 1, 0}; // advection coefficient
+        Eigen::Vector3d a = {1, 0, 0}; // advection coefficient
         Eigen::SparseMatrix<double> K;
         mesh.getStiffMatrix(K, a);
 
@@ -67,7 +67,7 @@ namespace solver {
         double t0 = 0;
         double tf = 1;
         int step = 0;
-        double dt = 0.1;
+        double dt = 0.01;
         for(double t=t0; t<tf; t+=dt, ++step) {
 
             // Savings