From 615d339f1cef7ea3f6d7cc9e62b5735afba432ea Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Fri, 9 Jul 2021 12:08:24 +0200
Subject: [PATCH] format
---
heat/src/wSolver.cpp | 297 ++++++++++++++++-----------------
mirrors/src/wSolver.cpp | 51 +++---
tbox/src/wMshDeform.cpp | 4 +-
waves/src/wTimeIntegration.cpp | 94 +++++------
4 files changed, 216 insertions(+), 230 deletions(-)
diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index a1890470..22a68553 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -256,48 +256,47 @@ 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";
-
- // (flux moyen . volume) sur l'element
- Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, *mat->k);
-
- 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';
- }
-
- tbb::spin_mutex::scoped_lock lock(mutex);
- Vtot += V;
- qM += qV; // qM = qM + qV
- });
+ 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);
+
+ 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';
+ }
+
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ Vtot += V;
+ qM += qV; // qM = qM + qV
+ });
}
qM /= Vtot;
@@ -372,26 +371,25 @@ 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";
-
- 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)
- {
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ke(ii, jj)));
- }
- }
- });
+ 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);
+
+ // 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>(rows[nodi->row], nodj->row, Ke(ii, jj)));
+ }
+ }
+ });
}
}
else
@@ -432,24 +430,23 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
if (verbose > 1)
std::cout << "\tprocessing " << *mat << '\n';
- std::for_each(mat->tag->elems.begin(), mat->tag->elems.end(), // devrait pouvoir etre fait en TBB...
- [&](Element *e) {
- if (e->type() != ELTYPE::TRI3)
- return;
- 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)
- {
- Node *nodj = e->nodes[jj];
- T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ke(ii, jj)));
- }
- }
- });
+ std::for_each(mat->tag->elems.begin(), mat->tag->elems.end(), [&](Element *e) { // devrait pouvoir etre fait en TBB...
+ if (e->type() != ELTYPE::TRI3)
+ return;
+ 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)
+ {
+ Node *nodj = e->nodes[jj];
+ T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ke(ii, jj)));
+ }
+ }
+ });
}
if (verbose > 1)
std::cout << "assembly done.\n";
@@ -524,26 +521,25 @@ 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";
-
- // ** 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);
-
- Vtot += Ce(jj);
- rhocM_Vtot += (mat->rhoc) * Ce(jj);
- }
- });
+ 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));
+
+ // 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);
+ }
+ });
}
rhs(n1->row) = rhocM_Vtot * TM;
if (verbose > 1)
@@ -596,22 +592,21 @@ 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";
-
- 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);
- }
- });
+ 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);
+
+ // 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())
@@ -646,24 +641,23 @@ 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";
-
- // ** se vector => s vector
- Eigen::VectorXd se = HeatTerm::build(*e, T1, *src->f);
-
- // 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);
- }
- });
+ 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);
+
+ // 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);
+ }
+ });
}
if (verbose > 1)
@@ -686,24 +680,23 @@ 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";
-
- // ** qe vector => q vector
- Eigen::VectorXd qe = HeatTerm::build(*e, T1, *bnd->f);
-
- // 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);
- }
- });
+ 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);
+
+ // 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);
+ }
+ });
}
if (verbose > 1)
diff --git a/mirrors/src/wSolver.cpp b/mirrors/src/wSolver.cpp
index 55f06b66..24ee49ed 100644
--- a/mirrors/src/wSolver.cpp
+++ b/mirrors/src/wSolver.cpp
@@ -163,33 +163,30 @@ void Solver::start(MshExport *mshWriter)
// K matrix
//=============================tbb==================================
- tbb::parallel_for_each(current_tag.elems.begin(), current_tag.elems.end(), [&](Element *e)
- //tbb::parallel_for_each(size_t(0), size_t(current_tag.elems.size()), [&](size_t index)
- //tbb::parallel_for_each(current_tag.elems.begin(), current_tag.elems.end(),[&](Element *e)
- {
- // Element *e = current_tag.elems[index];
-
- if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
- return; //if the element is not a tetrahedron or a hexahedron, it is not treated
-
- Eigen::MatrixXd K_e = ThermoMecaTerm::buildK(*e, H);
-
- 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];
- //=====================================//
- tbb::spin_mutex::scoped_lock lock(mutex);
- //=====================================//
- for (auto iii = 0; iii < 3; ++iii)
- for (auto jjj = 0; jjj < 3; ++jjj)
- if (K_e(3 * ii + iii, 3 * jj + jjj) != 0)
- Trip.push_back(Eigen::Triplet<double>(3 * (nodi->row) + iii, 3 * (nodj->row) + jjj, K_e(3 * ii + iii, 3 * jj + jjj)));
- }
- }
- });
+ tbb::parallel_for_each(current_tag.elems.begin(), current_tag.elems.end(), [&](Element *e) {
+ // Element *e = current_tag.elems[index];
+
+ if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
+ return; //if the element is not a tetrahedron or a hexahedron, it is not treated
+
+ Eigen::MatrixXd K_e = ThermoMecaTerm::buildK(*e, H);
+
+ 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];
+ //=====================================//
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ //=====================================//
+ for (auto iii = 0; iii < 3; ++iii)
+ for (auto jjj = 0; jjj < 3; ++jjj)
+ if (K_e(3 * ii + iii, 3 * jj + jjj) != 0)
+ Trip.push_back(Eigen::Triplet<double>(3 * (nodi->row) + iii, 3 * (nodj->row) + jjj, K_e(3 * ii + iii, 3 * jj + jjj)));
+ }
+ }
+ });
//=============================tbb==================================
}
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index f10b2ccb..22463cc4 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -302,7 +302,7 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
std::deque<Eigen::Triplet<double>> T;
// Build stiffness matrix
- tbb::parallel_for_each(fldElems.begin(), fldElems.end(), [&](Element* e) {
+ tbb::parallel_for_each(fldElems.begin(), fldElems.end(), [&](Element *e) {
// Hooke tensor
Eigen::MatrixXd H;
if (nDim == 2)
@@ -336,10 +336,8 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
H(2, 1) = H(1, 2);
H *= E / ((1 + nu) * (1 - 2 * nu));
}
-
// Elementary stiffness matrix
Eigen::MatrixXd Ke = this->buildK(*e, H);
-
// Assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < e->nodes.size(); ++i)
diff --git a/waves/src/wTimeIntegration.cpp b/waves/src/wTimeIntegration.cpp
index f96c7a3b..d5847512 100644
--- a/waves/src/wTimeIntegration.cpp
+++ b/waves/src/wTimeIntegration.cpp
@@ -102,28 +102,27 @@ 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";
-
- // ** Se matrix => S vector
- Eigen::MatrixXd Se = WaveTerm::buildM(*e);
-
- // 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, bnd->c * Se(ii, jj)));
- }
- }
- });
+ 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);
+
+ // 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, bnd->c * Se(ii, jj)));
+ }
+ }
+ });
}
// Build, clean and turn to compressed row format
S2.setFromTriplets(T.begin(), T.end());
@@ -159,32 +158,31 @@ 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";
-
- // ** Me matrix => Md vector
- Eigen::MatrixXd Me = WaveTerm::buildM(*e);
-
- // ** Ke matrix => K matrix
- Eigen::MatrixXd Ke = WaveTerm::buildK(*e, u);
-
- // 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, (mat->c * mat->c) * Ke(ii, jj)));
- Md[nodi->row] += Me(ii, jj);
- }
- }
- });
+ 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);
+
+ // ** Ke matrix => K matrix
+ Eigen::MatrixXd Ke = WaveTerm::buildK(*e, u);
+
+ // 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, (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());
--
GitLab