diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index 96c7f606d8ac643c99e2403ff2751f8d8c39bd9e..30e65f0e610d52b14f6e32d85d26e997b904fbab 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -256,47 +256,50 @@ void Solver::start(MshExport *mshWriter)
Eigen::Vector2d qM = Eigen::Vector2d::Zero();
for (auto mat : pbl->media)
{
- tbb::parallel_for_each(mat->tag->elems.begin(), mat->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::TRI3)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ mat->tag->elems.begin(), mat->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::TRI3)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- // (flux moyen . volume) sur l'element
- Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, *mat->k);
+ // (flux moyen . volume) sur l'element
+ Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, *mat->k);
- double V = e->getVol();
+ double V = e->getVol();
- if (save)
- {
- // gradient of T at GP 0
- Eigen::VectorXd grad = e->computeGradient(T1, 0);
- int i = e->no - 1;
- gradT_x[i] = grad(0);
- gradT_y[i] = grad(1);
- gradT[i] = Eigen::Vector3d(grad(0), grad(1), 0);
-
- // mean q over the element
- kgradT[i] = Eigen::Vector3d(qV(0) / V, qV(1) / V, 0);
-
- // mean k over the element
- Eigen::MatrixXd kmean = HeatTerm::computeMatrix(*e, T1, *mat->k) / V;
- kmoy11[i] = kmean(0, 0);
- kmoy22[i] = kmean(1, 1);
- kmoy12[i] = kmean(0, 1);
- kmoy21[i] = kmean(1, 0);
-
- // tensor
- kmean.conservativeResize(3, 3);
- kmean.row(2) = Eigen::RowVector3d::Zero();
- kmean.col(2) = Eigen::Vector3d::Zero();
- kmoy[i] = kmean;
- //std::cout << kmean << '\n';
- }
+ if (save)
+ {
+ // gradient of T at GP 0
+ Eigen::VectorXd grad = e->computeGradient(T1, 0);
+ int i = e->no - 1;
+ gradT_x[i] = grad(0);
+ gradT_y[i] = grad(1);
+ gradT[i] = Eigen::Vector3d(grad(0), grad(1), 0);
+
+ // mean q over the element
+ kgradT[i] = Eigen::Vector3d(qV(0) / V, qV(1) / V, 0);
+
+ // mean k over the element
+ Eigen::MatrixXd kmean = HeatTerm::computeMatrix(*e, T1, *mat->k) / V;
+ kmoy11[i] = kmean(0, 0);
+ kmoy22[i] = kmean(1, 1);
+ kmoy12[i] = kmean(0, 1);
+ kmoy21[i] = kmean(1, 0);
+
+ // tensor
+ kmean.conservativeResize(3, 3);
+ kmean.row(2) = Eigen::RowVector3d::Zero();
+ kmean.col(2) = Eigen::Vector3d::Zero();
+ kmoy[i] = kmean;
+ // std::cout << kmean << '\n';
+ }
- tbb::spin_mutex::scoped_lock lock(mutex);
- Vtot += V;
- qM += qV; // qM = qM + qV });
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ Vtot += V;
+ qM += qV; // qM = qM + qV
+ });
}
qM /= Vtot;
@@ -371,25 +374,28 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
{
if (verbose > 1)
std::cout << "\tprocessing " << *mat << '\n';
- tbb::parallel_for_each(mat->tag->elems.begin(), mat->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::TRI3)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ mat->tag->elems.begin(), mat->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::TRI3)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, *mat->k, false);
+ Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, *mat->k, false);
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ke(ii, jj)));
+ Node *nodi = e->nodes[ii];
+ for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ {
+ Node *nodj = e->nodes[jj];
+ T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ke(ii, jj)));
+ }
}
- } });
+ });
}
}
else
@@ -522,25 +528,28 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
for (auto mat : pbl->media)
{
// std::cout << "\tprocessing " << *mat << '\n';
- tbb::parallel_for_each(mat->tag->elems.begin(), mat->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::TRI3)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ mat->tag->elems.begin(), mat->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::TRI3)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- // ** Ce matrix => K matrix
- Eigen::VectorXd Ce = HeatTerm::build(*e, T1, Fct0C(1.0));
+ // ** Ce matrix => K matrix
+ Eigen::VectorXd Ce = HeatTerm::build(*e, T1, Fct0C(1.0));
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t jj = 0; jj < e->nodes.size(); ++jj)
- {
- Node *nodj = e->nodes[jj];
- K2.coeffRef(n1->row, nodj->row) += (mat->rhoc) * Ce(jj);
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ {
+ Node *nodj = e->nodes[jj];
+ K2.coeffRef(n1->row, nodj->row) += (mat->rhoc) * Ce(jj);
- Vtot += Ce(jj);
- rhocM_Vtot += (mat->rhoc) * Ce(jj);
- } });
+ Vtot += Ce(jj);
+ rhocM_Vtot += (mat->rhoc) * Ce(jj);
+ }
+ });
}
rhs(n1->row) = rhocM_Vtot * TM;
if (verbose > 1)
@@ -593,21 +602,24 @@ void Solver::buildqint(Eigen::VectorXd &qint)
{
if (verbose > 1)
std::cout << "\tprocessing " << *mat << '\n';
- tbb::parallel_for_each(mat->tag->elems.begin(), mat->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::TRI3)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ mat->tag->elems.begin(), mat->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::TRI3)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- Eigen::VectorXd qe = HeatTerm::build2(*e, T1, *mat->k);
+ Eigen::VectorXd qe = HeatTerm::build2(*e, T1, *mat->k);
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- qint(nodi->row) += qe(ii);
- } });
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
+ {
+ Node *nodi = e->nodes[ii];
+ qint(nodi->row) += qe(ii);
+ }
+ });
}
if (!pbl->dBCs.empty())
@@ -642,23 +654,26 @@ void Solver::builds(Eigen::Map<Eigen::VectorXd> &s)
for (auto src : pbl->srcs)
{
// std::cout << "\tprocessing " << *src << '\n';
- tbb::parallel_for_each(src->tag->elems.begin(), src->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::TRI3 && e->type() != ElType::LINE2)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ src->tag->elems.begin(), src->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::TRI3 && e->type() != ElType::LINE2)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- // ** se vector => s vector
- Eigen::VectorXd se = HeatTerm::build(*e, T1, *src->f);
+ // ** se vector => s vector
+ Eigen::VectorXd se = HeatTerm::build(*e, T1, *src->f);
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- s(nodi->row) += se(ii);
- } });
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
+ {
+ Node *nodi = e->nodes[ii];
+ s(nodi->row) += se(ii);
+ }
+ });
}
if (verbose > 1)
@@ -681,23 +696,26 @@ void Solver::buildq(Eigen::Map<Eigen::VectorXd> &q)
for (auto bnd : pbl->bnds)
{
// std::cout << "\tprocessing " << *src << '\n';
- tbb::parallel_for_each(bnd->tag->elems.begin(), bnd->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::LINE2)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ bnd->tag->elems.begin(), bnd->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::LINE2)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- // ** qe vector => q vector
- Eigen::VectorXd qe = HeatTerm::build(*e, T1, *bnd->f);
+ // ** qe vector => q vector
+ Eigen::VectorXd qe = HeatTerm::build(*e, T1, *bnd->f);
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- q(nodi->row) += qe(ii);
- } });
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
+ {
+ Node *nodi = e->nodes[ii];
+ q(nodi->row) += qe(ii);
+ }
+ });
}
if (verbose > 1)
diff --git a/waves/src/wTimeIntegration.cpp b/waves/src/wTimeIntegration.cpp
index 074f7bd9a440a38dc2791504263cdcc373f801af..83f5aa86c67e056b60583ad9c4e60cfebeacfd87 100644
--- a/waves/src/wTimeIntegration.cpp
+++ b/waves/src/wTimeIntegration.cpp
@@ -102,27 +102,30 @@ void TimeIntegration::buildS(Eigen::SparseMatrix<double, Eigen::RowMajor> &S2)
{
std::cout << "\tprocessing " << *bnd << '\n';
std::cout << "normal=" << static_cast<Quad4 *>(bnd->tag->elems[0])->getNormal() << '\n';
- tbb::parallel_for_each(bnd->tag->elems.begin(), bnd->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::QUAD4)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ bnd->tag->elems.begin(), bnd->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::QUAD4)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- // ** Se matrix => S vector
- Eigen::MatrixXd Se = WaveTerm::buildM(*e);
+ // ** Se matrix => S vector
+ Eigen::MatrixXd Se = WaveTerm::buildM(*e);
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, bnd->c * Se(ii, jj)));
+ Node *nodi = e->nodes[ii];
+ for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ {
+ Node *nodj = e->nodes[jj];
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, bnd->c * Se(ii, jj)));
+ }
}
- } });
+ });
}
// Build, clean and turn to compressed row format
S2.setFromTriplets(T.begin(), T.end());
@@ -158,31 +161,34 @@ void TimeIntegration::buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowM
for (auto mat : pbl->media)
{
std::cout << "\tprocessing " << *mat << '\n';
- tbb::parallel_for_each(mat->tag->elems.begin(), mat->tag->elems.end(), [&](Element *e)
- {
- if (e->type() != ElType::HEX8)
- return;
- //std::cout << "processing element #" << e->no << "\n";
+ tbb::parallel_for_each(
+ mat->tag->elems.begin(), mat->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::HEX8)
+ return;
+ // std::cout << "processing element #" << e->no << "\n";
- // ** Me matrix => Md vector
- Eigen::MatrixXd Me = WaveTerm::buildM(*e);
+ // ** Me matrix => Md vector
+ Eigen::MatrixXd Me = WaveTerm::buildM(*e);
- // ** Ke matrix => K matrix
- Eigen::MatrixXd Ke = WaveTerm::buildK(*e, u);
+ // ** Ke matrix => K matrix
+ Eigen::MatrixXd Ke = WaveTerm::buildK(*e, u);
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, (mat->c * mat->c) * Ke(ii, jj)));
- Md[nodi->row] += Me(ii, jj);
+ Node *nodi = e->nodes[ii];
+ for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ {
+ Node *nodj = e->nodes[jj];
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, (mat->c * mat->c) * Ke(ii, jj)));
+ Md[nodi->row] += Me(ii, jj);
+ }
}
- } });
+ });
}
// Build, clean and turn to compressed row format
K2.setFromTriplets(T.begin(), T.end());
@@ -194,7 +200,8 @@ void TimeIntegration::buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowM
std::cout << "[cpu] " << chrono1 << '\n';
}
-/*void TimeIntegration::build(MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2)
+/*
+void TimeIntegration::build(MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2)
{
tbb::spin_mutex mutex;
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads); // TODO mettre ça ailleurs...
@@ -207,24 +214,26 @@ void TimeIntegration::buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowM
for (auto mat : pbl->media)
{
std::cout << "\tprocessing " << *mat << '\n';
- tbb::parallel_for_each(mat->tag->elems.begin(), mat->tag->elems.end(),
- [&](Element *e) {
- if (e->type() != ElType::HEX8)
- return;
- Eigen::MatrixXd Ae = e->build(type);
-
- // assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t ii = 0; ii < e->nodes.size(); ++ii)
- {
- Node *nodi = e->nodes[ii];
- for (size_t jj = 0; jj < e->nodes.size(); ++jj)
- {
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae(ii, jj)));
- }
- }
- });
+ tbb::parallel_for_each(
+ mat->tag->elems.begin(), mat->tag->elems.end(),
+ [&](Element *e)
+ {
+ if (e->type() != ElType::HEX8)
+ return;
+ Eigen::MatrixXd Ae = e->build(type);
+
+ // assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t ii = 0; ii < e->nodes.size(); ++ii)
+ {
+ Node *nodi = e->nodes[ii];
+ for (size_t jj = 0; jj < e->nodes.size(); ++jj)
+ {
+ Node *nodj = e->nodes[jj];
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae(ii, jj)));
+ }
+ }
+ });
}
// Build, clean and turn to compressed row format
A2.setFromTriplets(T.begin(), T.end());
@@ -234,7 +243,8 @@ void TimeIntegration::buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowM
std::cout << "S (" << A2.rows() << "," << A2.cols() << ") nnz=" << A2.nonZeros() << "\n";
chrono1.read();
std::cout << "[cpu] " << chrono1 << '\n';
-}*/
+}
+*/
void TimeIntegration::write(std::ostream &out) const
{