diff --git a/srcs/BEM/mainBEM.cpp b/srcs/BEM/mainBEM.cpp index 02035781a5f8ac07e30bc7e0ac8848a8994c4116..1174d1849e89f97d56c82b214f4dba066fd5f90f 100644 --- a/srcs/BEM/mainBEM.cpp +++ b/srcs/BEM/mainBEM.cpp @@ -11,7 +11,10 @@ int main(int argc, char **argv) return 0; } + double start = omp_get_wtime(); + int nthreads = omp_get_max_threads(); gmsh::initialize(); + omp_set_num_threads(nthreads); gmsh::open(argv[1]); gmsh::model::mesh::generate(2); diff --git a/srcs/BEM/postProcessing.cpp b/srcs/BEM/postProcessing.cpp index c5854559391eff43301a958a05c5861d9c752c03..408d57b0bd40473218dcbc3576f4fcf8a900c28b 100644 --- a/srcs/BEM/postProcessing.cpp +++ b/srcs/BEM/postProcessing.cpp @@ -377,6 +377,7 @@ void computeElementData(std::vector<std::vector<std::vector<std::vector<double>> phi[i][j].resize(domainElementTags[i][j].size()); electricField[i][j].resize(domainElementTags[i][j].size()); gmsh::model::mesh::getBarycenters(domainElementTypes[i][j], domainEntityTags[i], false, false, barycenters); + #pragma omp parallel for schedule(static) for(size_t k = 0; k < domainElementTags[i][j].size(); k++) { double x = barycenters[3*k]; @@ -386,55 +387,44 @@ void computeElementData(std::vector<std::vector<std::vector<std::vector<double>> double gradPhiX = 0; double gradPhiY = 0; - #pragma omp parallel shared(currentPhi, gradPhiX, gradPhiY) + double localCurrentPhi = 0; + double localGradPhiX = 0; + double localGradPhiY = 0; + + double h, hX, hY; + double g, gX, gY; + + for(std::size_t el = 0; el < elementVector.size(); el++) { - double localCurrentPhi = 0; - double localGradPhiX = 0; - double localGradPhiY = 0; - - double h, hX, hY; - double g, gX, gY; - - #pragma omp for schedule(guided) - for(std::size_t el = 0; el < elementVector.size(); el++) + if(computePhi) { - if(computePhi) - { - h = computeAnalyticalH(x, y, elementVector[el]); - if(!solutionType) - g = computeNumericalG(x, y, elementVector[el], integrationType); - else - g = computeAnalyticalG(x, y, elementVector[el]); + h = computeAnalyticalH(x, y, elementVector[el]); + if(!solutionType) + g = computeNumericalG(x, y, elementVector[el], integrationType); + else + g = computeAnalyticalG(x, y, elementVector[el]); - localCurrentPhi += h*elementVector[el].bcValue[0] - g*elementVector[el].bcValue[1]; + localCurrentPhi += h*elementVector[el].bcValue[0] - g*elementVector[el].bcValue[1]; + } + if(computeElectricField) + { + if(!solutionType) + { + hX = computeNumericalGradHx(x, y, elementVector[el], integrationType); + gX = computeNumericalGradGx(x, y, elementVector[el], integrationType); + hY = computeNumericalGradHy(x, y, elementVector[el], integrationType); + gY = computeNumericalGradGy(x, y, elementVector[el], integrationType); } - if(computeElectricField) + else { - if(!solutionType) - { - hX = computeNumericalGradHx(x, y, elementVector[el], integrationType); - gX = computeNumericalGradGx(x, y, elementVector[el], integrationType); - hY = computeNumericalGradHy(x, y, elementVector[el], integrationType); - gY = computeNumericalGradGy(x, y, elementVector[el], integrationType); - } - else - { - hX = computeAnalyticalGradHx(x, y, elementVector[el]); - gX = computeAnalyticalGradGx(x, y, elementVector[el]); - hY = computeAnalyticalGradHy(x, y, elementVector[el]); - gY = computeAnalyticalGradGy(x, y, elementVector[el]); - } - - localGradPhiX += hX*elementVector[el].bcValue[0] - gX*elementVector[el].bcValue[1]; - localGradPhiY += hY*elementVector[el].bcValue[0] - gY*elementVector[el].bcValue[1]; + hX = computeAnalyticalGradHx(x, y, elementVector[el]); + gX = computeAnalyticalGradGx(x, y, elementVector[el]); + hY = computeAnalyticalGradHy(x, y, elementVector[el]); + gY = computeAnalyticalGradGy(x, y, elementVector[el]); } - } - #pragma omp critical - { - currentPhi += localCurrentPhi; - gradPhiX += localGradPhiX; - gradPhiY += localGradPhiY; + localGradPhiX += hX*elementVector[el].bcValue[0] - gX*elementVector[el].bcValue[1]; + localGradPhiY += hY*elementVector[el].bcValue[0] - gY*elementVector[el].bcValue[1]; } }