Boman Romain
--- a/flow/src/wBoundary.cpp

+ 11

− 44
+++ b/flow/src/wBoundary.cpp

+ 11

− 44
 @@ -96,55 +96,22 @@ void Boundary::create()
        svMap[f->el0] = f->el1;
        delete f;
    }
+    // Resize load vectors
+    cLoadX.resize(nodes.size());
+    cLoadY.resize(nodes.size());
+    cLoadZ.resize(nodes.size());
 }

 /*
- * @brief Integrate f over elements and extract
+ * @brief Integrate f over an element
 * @authors Adrien Crovato
 */
-std::vector<double> Boundary::integrate(std::vector<double> const &u, Fct0 const &f) const
+double Boundary::integrate(Element *e, std::vector<double> const &u, tbox::Fct0 const &f) const
 {
-    std::vector<double> vI(svMap.size());
-    for (size_t i = 0; i < vI.size(); ++i)
-    {
-        Element *e = groups[0]->tag->elems[i];
-        // evaluate the function in the volume element
-        double fk = f.eval(*svMap.at(e), u, 0);
-        // integrate the function on the surface element
-        Gauss &gauss = e->getVCache().getVGauss();
-        for (size_t k = 0; k < gauss.p.size(); ++k)
-            vI[i] += fk * gauss.w[k] * e->getVMem().getDetJ(k);
-    }
-    return vI;
-}
-
-/**
- * @brief Transfer u from elements to nodes and extract (vector)
- * @authors Adrien Crovato
- */
-std::vector<std::vector<double>> Boundary::transfer(std::vector<double> const &u) const
-{
-    // Map each element to its value
-    std::map<Element *, double> euMap;
-    for (size_t i = 0; i < groups[0]->tag->elems.size(); i++)
-        euMap[groups[0]->tag->elems[i]] = u[i];
-
-    // Return a vector (x, y, z) of summed neighbour values for each node
-    // This is a "convervative" way of transfering u as the sum of u over the
-    // elements will the same as the sum of u over the nodes
-    std::vector<std::vector<double>> vN(nodes.size());
-    for (size_t i = 0; i < vN.size(); ++i)
-    {
-        std::vector<double> vInt(3, 0.);
-        for (auto e : neMap.at(nodes[i]))
-        {
-            vInt[0] += euMap.at(e) * e->normal().x[0] / e->nodes.size();
-            vInt[1] += euMap.at(e) * e->normal().x[1] / e->nodes.size();
-            vInt[2] += euMap.at(e) * e->normal().x[2] / e->nodes.size();
-        }
-        vN[i] = vInt;
-    }
-    return vN;
+    // evaluate the function in the volume element
+    double fk = f.eval(*svMap.at(e), u, 0);
+    // integrate the function on the surface element
+    return fk * e->getVMem().getVol(); // since fk is constant: sum(fk*wk*dtmk) = fk*V
 }

 /**
 @@ -171,7 +138,7 @@ void Boundary::save(std::string const &name, Results const &res)
            << std::setw(15) << Cl
            << std::setw(15) << Cd
            << std::setw(15) << Cm
-            << std::endl;       
+            << std::endl;
    // Elements (connectvity)
    outfile << "$Elements" << std::endl;
    outfile << groups[0]->tag->elems.size() << std::endl;