@@ -2338,7 +2338,7 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config)
23382338 /* --- Update the control volume structures ---*/
23392339
23402340 geometry_container[iMesh]->SetControlVolume (geometry_container[iMesh - 1 ], UPDATE);
2341- geometry_container[iMesh]->SetBoundControlVolume (geometry_container[iMesh - 1 ], UPDATE);
2341+ geometry_container[iMesh]->SetBoundControlVolume (geometry_container[iMesh - 1 ], config, UPDATE);
23422342 geometry_container[iMesh]->SetCoord (geometry_container[iMesh - 1 ]);
23432343 }
23442344
@@ -2453,123 +2453,71 @@ su2double CGeometry::GetSurfaceArea(const CConfig* config, unsigned short val_ma
24532453 return 0.0 ;
24542454}
24552455
2456- void CGeometry::ComputeSurf_Straightness (CConfig* config, bool print_on_screen) {
2457- bool RefUnitNormal_defined;
2458- unsigned short iDim, iMarker, iMarker_Global, nMarker_Global = config->GetnMarker_CfgFile ();
2459- unsigned long iVertex;
2460- constexpr passivedouble epsilon = 1.0e-6 ;
2461- su2double Area;
2462- string Local_TagBound, Global_TagBound;
2463-
2464- vector<su2double> Normal (nDim), UnitNormal (nDim), RefUnitNormal (nDim);
2465-
2466- /* --- Assume now that this boundary marker is straight. As soon as one
2467- AreaElement is found that is not aligend with a Reference then it is
2468- certain that the boundary marker is not straight and one can stop
2469- searching. Another possibility is that this process doesn't own
2470- any nodes of that boundary, in that case we also have to assume the
2471- boundary is straight.
2472- Any boundary type other than SYMMETRY_PLANE or EULER_WALL gets
2473- the value false (or see cases specified in the conditional below)
2474- which could be wrong. ---*/
2475- bound_is_straight.resize (nMarker);
2476- fill (bound_is_straight.begin (), bound_is_straight.end (), true );
2477-
2478- /* --- Loop over all local markers ---*/
2479- for (iMarker = 0 ; iMarker < nMarker; iMarker++) {
2480- Local_TagBound = config->GetMarker_All_TagBound (iMarker);
2481-
2482- /* --- Marker has to be Symmetry or Euler. Additionally marker can't be a
2483- moving surface and Grid Movement Elasticity is forbidden as well. All
2484- other GridMovements are rigid. ---*/
2485- if ((config->GetMarker_All_KindBC (iMarker) == SYMMETRY_PLANE ||
2486- config->GetMarker_All_KindBC (iMarker) == EULER_WALL) &&
2487- !config->GetMarker_Moving_Bool (Local_TagBound) && !config->GetMarker_Deform_Mesh_Bool (Local_TagBound)) {
2488- /* --- Loop over all global markers, and find the local-global pair via
2489- matching unique string tags. ---*/
2490- for (iMarker_Global = 0 ; iMarker_Global < nMarker_Global; iMarker_Global++) {
2491- Global_TagBound = config->GetMarker_CfgFile_TagBound (iMarker_Global);
2492- if (Local_TagBound == Global_TagBound) {
2493- RefUnitNormal_defined = false ;
2494- iVertex = 0 ;
2495-
2496- while (bound_is_straight[iMarker] && iVertex < nVertex[iMarker]) {
2497- vertex[iMarker][iVertex]->GetNormal (Normal.data ());
2498- UnitNormal = Normal;
2499-
2500- /* --- Compute unit normal. ---*/
2501- Area = 0.0 ;
2502- for (iDim = 0 ; iDim < nDim; iDim++) Area += Normal[iDim] * Normal[iDim];
2503- Area = sqrt (Area);
2504-
2505- /* --- Negate for outward convention. ---*/
2506- for (iDim = 0 ; iDim < nDim; iDim++) UnitNormal[iDim] /= -Area;
2507-
2508- /* --- Check if unit normal is within tolerance of the Reference unit normal.
2509- Reference unit normal = first unit normal found. ---*/
2510- if (RefUnitNormal_defined) {
2511- for (iDim = 0 ; iDim < nDim; iDim++) {
2512- if (abs (RefUnitNormal[iDim] - UnitNormal[iDim]) > epsilon) {
2513- bound_is_straight[iMarker] = false ;
2514- break ;
2515- }
2516- }
2517- } else {
2518- RefUnitNormal = UnitNormal; // deep copy of values
2519- RefUnitNormal_defined = true ;
2520- }
2456+ void CGeometry::ComputeModifiedSymmetryNormals (const CConfig* config) {
2457+ /* Check how many symmetry planes there are and use the first (lowest ID) as the basis to orthogonalize against.
2458+ * All nodes that are shared by multiple symmetries have to get a corrected normal. */
2459+ symmetryNormals.resize (nMarker);
2460+ std::vector<unsigned short > symMarkers;
25212461
2522- iVertex++;
2523- } // while iVertex
2524- } // if Local == Global
2525- } // for iMarker_Global
2526- } else {
2527- /* --- Enforce default value: false ---*/
2528- bound_is_straight[iMarker] = false ;
2529- } // if sym or euler ...
2530- } // for iMarker
2531-
2532- /* --- Communicate results and print on screen. ---*/
2533- if (print_on_screen) {
2534- /* --- Additional vector which can later be MPI::Allreduce(d) to pring the results
2535- on screen as nMarker (local) can vary across ranks. Default 'true' as it can
2536- happen that a local rank does not contain an element of each surface marker. ---*/
2537- vector<bool > bound_is_straight_Global (nMarker_Global, true );
2538- /* --- Match local with global tag bound and fill a Global Marker vector. ---*/
2539- for (iMarker = 0 ; iMarker < nMarker; iMarker++) {
2540- Local_TagBound = config->GetMarker_All_TagBound (iMarker);
2541- for (iMarker_Global = 0 ; iMarker_Global < nMarker_Global; iMarker_Global++) {
2542- Global_TagBound = config->GetMarker_CfgFile_TagBound (iMarker_Global);
2543-
2544- if (Local_TagBound == Global_TagBound) bound_is_straight_Global[iMarker_Global] = bound_is_straight[iMarker];
2545-
2546- } // for iMarker_Global
2547- } // for iMarker
2548-
2549- vector<int > Buff_Send_isStraight (nMarker_Global), Buff_Recv_isStraight (nMarker_Global);
2550-
2551- /* --- Cast to int as std::vector<boolean> can be a special construct. MPI handling using <int>
2552- is more straight-forward. ---*/
2553- for (iMarker_Global = 0 ; iMarker_Global < nMarker_Global; iMarker_Global++)
2554- Buff_Send_isStraight[iMarker_Global] = static_cast <int >(bound_is_straight_Global[iMarker_Global]);
2555-
2556- /* --- Product of type <int>(bool) is equivalnt to a 'logical and' ---*/
2557- SU2_MPI::Allreduce (Buff_Send_isStraight.data (), Buff_Recv_isStraight.data (), nMarker_Global, MPI_INT, MPI_PROD,
2558- SU2_MPI::GetComm ());
2559-
2560- /* --- Print results on screen. ---*/
2561- if (rank == MASTER_NODE) {
2562- for (iMarker_Global = 0 ; iMarker_Global < nMarker_Global; iMarker_Global++) {
2563- if (config->GetMarker_CfgFile_KindBC (config->GetMarker_CfgFile_TagBound (iMarker_Global)) == SYMMETRY_PLANE ||
2564- config->GetMarker_CfgFile_KindBC (config->GetMarker_CfgFile_TagBound (iMarker_Global)) == EULER_WALL) {
2565- cout << " Boundary marker " GetMarker_CfgFile_TagBound (iMarker_Global) << " is"
2566- if (!static_cast <bool >(Buff_Recv_isStraight[iMarker_Global])) cout << " NOT"
2567- if (nDim == 2 ) cout << " a single straight."
2568- if (nDim == 3 ) cout << " a single plane."
2569- } // if sym or euler
2570- } // for iMarker_Global
2571- } // if rank==MASTER
2572- } // if print_on_scren
2462+ for (auto iMarker = 0u ; iMarker < nMarker; ++iMarker) {
2463+ if ((config->GetMarker_All_KindBC (iMarker) == SYMMETRY_PLANE) ||
2464+ (config->GetMarker_All_KindBC (iMarker) == EULER_WALL)) {
2465+ symMarkers.push_back (iMarker);
2466+ }
2467+ }
2468+
2469+ /* --- Loop over all markers and find nodes on symmetry planes that are shared with other symmetries. ---*/
2470+ /* --- The first symmetry does not need a corrected normal vector, hence start at 1. ---*/
2471+ for (size_t i = 1 ; i < symMarkers.size (); ++i) {
2472+ const auto iMarker = symMarkers[i];
2473+
2474+ for (auto iVertex = 0ul ; iVertex < nVertex[iMarker]; iVertex++) {
2475+ const auto iPoint = vertex[iMarker][iVertex]->GetNode ();
2476+
2477+ /* --- Halo points do not need to be considered. ---*/
2478+ if (!nodes->GetDomain (iPoint)) continue ;
2479+
2480+ /* --- Get the vertex normal on the current symmetry. ---*/
2481+ std::array<su2double, MAXNDIM> iNormal = {};
2482+ vertex[iMarker][iVertex]->GetNormal (iNormal.data ());
2483+
2484+ /* --- Loop over previous symmetries and if this point shares them, make this normal orthogonal to them. ---*/
2485+ bool isShared = false ;
2486+
2487+ for (size_t j = 0 ; j < i; ++j) {
2488+ const auto jMarker = symMarkers[j];
2489+ const auto jVertex = nodes->GetVertex (iPoint, jMarker);
2490+ if (jVertex < 0 ) continue ;
2491+
2492+ isShared = true ;
2493+
2494+ std::array<su2double, MAXNDIM> jNormal = {};
2495+ const auto it = symmetryNormals[jMarker].find (jVertex);
2496+
2497+ if (it != symmetryNormals[jMarker].end ()) {
2498+ jNormal = it->second ;
2499+ } else {
2500+ vertex[jMarker][jVertex]->GetNormal (jNormal.data ());
2501+ const su2double area = GeometryToolbox::Norm (nDim, jNormal.data ());
2502+ for (auto iDim = 0u ; iDim < nDim; iDim++) jNormal[iDim] /= area;
2503+ }
2504+
2505+ const auto proj = GeometryToolbox::DotProduct (nDim, jNormal.data (), iNormal.data ());
2506+ for (auto iDim = 0u ; iDim < nDim; iDim++) iNormal[iDim] -= proj * jNormal[iDim];
2507+ }
2508+
2509+ if (!isShared) continue ;
2510+
2511+ /* --- Normalize. If the norm is close to zero it means the normal is a linear combination of previous
2512+ * normals, in this case we don't need to store the corrected normal, using the original in the gradient
2513+ * correction will have no effect since previous corrections will remove components in this direction). ---*/
2514+ const su2double area = GeometryToolbox::Norm (nDim, iNormal.data ());
2515+ if (area > EPS) {
2516+ for (auto iDim = 0u ; iDim < nDim; iDim++) iNormal[iDim] /= area;
2517+ symmetryNormals[iMarker][iVertex] = iNormal;
2518+ }
2519+ }
2520+ }
25732521}
25742522
25752523void CGeometry::ComputeSurf_Curvature (CConfig* config) {
0 commit comments