1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
|
// Copyright (C) 2025 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#include <iostream>
#include <algorithm>
#include <Eigen/Dense>
#include "StOpt/core/grids/GridAdapt2D.h"
#include "StOpt/core/utils/constant.h"
using namespace Eigen;
using namespace std;
using namespace StOpt;
GridAdapt2D::GridAdapt2D(const double &p_xMin, const double &p_xMax, const int &p_nbMeshX, const double &p_yMin, const double &p_yMax, const int &p_nbMeshY):
m_xMin(p_xMin), m_xMax(p_xMax), m_yMin(p_yMin), m_yMax(p_yMax)
{
double dx = (p_xMax - p_xMin) / p_nbMeshX;
double dy = (p_yMax - p_yMin) / p_nbMeshY;
for (int j = 0; j < p_nbMeshY; ++j)
{
double yB = p_yMin + j * dy;
for (int i = 0; i < p_nbMeshX ; ++i)
{
double xL = p_xMin + i * dx;
shared_ptr<Mesh2D> theMesh = make_shared<Mesh2D>(xL, xL + dx, yB, yB + dy, (j + 1) * (p_nbMeshX + 1) + i, j * (p_nbMeshX + 1) + i, j * (p_nbMeshX + 1) + i + 1, (j + 1) * (p_nbMeshX + 1) + i + 1);
vector< ArrayXi > vecPos = {{{i, j}}} ; //{aPos};
shared_ptr< vector< ArrayXi > > thePosition = make_shared<vector< ArrayXi > >(vecPos);
m_meshes.push_back(make_pair(theMesh, thePosition));
}
}
// store points
m_points.reserve((p_nbMeshX + 1) * (p_nbMeshY + 1));
for (int j = 0; j < p_nbMeshY + 1; ++j)
{
double y = p_yMin + j * dy;
for (int i = 0; i < p_nbMeshX + 1 ; ++i)
{
double x = p_xMin + i * dx;
m_points.push_back({{x, y}});
}
}
}
GridAdapt2D::GridAdapt2D(const double &p_xMin, const double &p_xMax, const double &p_yMin, const double &p_yMax, \
const list< pair< shared_ptr< Mesh2D>, std::shared_ptr< vector< ArrayXi > > > > &p_meshes, const vector< ArrayXd > &p_points) : GridAdaptBase(p_points), m_meshes(p_meshes), m_xMin(p_xMin), m_xMax(p_xMax), m_yMin(p_yMin), m_yMax(p_yMax)
{}
GridAdapt2D::GridAdapt2D(const double &p_xMin, const double &p_xMax, const double &p_yMin, const double &p_yMax, \
const vector< ArrayXd > &p_points) : GridAdaptBase(p_points), m_xMin(p_xMin), m_xMax(p_xMax), m_yMin(p_yMin), m_yMax(p_yMax)
{
}
void GridAdapt2D::splitMesh(pair< shared_ptr< Mesh2D >, std::shared_ptr<vector<ArrayXi > > > &p_meshToSplit)
{
// current number of points on grid
size_t nbpt = m_points.size();
pair<array< shared_ptr< Mesh2D >, 4 >, vector<ArrayXd > > newMeshAndPoints = p_meshToSplit.first->split(m_points);
std::shared_ptr<vector<ArrayXi > > iposToAdapt = p_meshToSplit.second;
// erase old mesh
auto it = find(m_meshes.begin(), m_meshes.end(), p_meshToSplit);
m_meshes.erase(it);
// store new meshes
// mesh splitting
// 0 3
// 1 2
std::shared_ptr<vector<ArrayXi > > iposToAdapt1 = make_shared<vector<ArrayXi > >(*iposToAdapt);
iposToAdapt1->push_back({{0, 1}}); // (0,1) position in square
m_meshes.push_back(make_pair(newMeshAndPoints.first[0], iposToAdapt1));
std::shared_ptr<vector<ArrayXi > > iposToAdapt2 = make_shared<vector<ArrayXi > >(*iposToAdapt);
iposToAdapt2->push_back({{0, 0}}); // (0,0) position in square
m_meshes.push_back(make_pair(newMeshAndPoints.first[1], iposToAdapt2));
std::shared_ptr<vector<ArrayXi > > iposToAdapt3 = make_shared<vector<ArrayXi > >(*iposToAdapt);
iposToAdapt3->push_back({{1, 0}}); // (1,0) position in square
m_meshes.push_back(make_pair(newMeshAndPoints.first[2], iposToAdapt3));
std::shared_ptr<vector<ArrayXi > > iposToAdapt4 = make_shared<vector<ArrayXi > >(*iposToAdapt);
iposToAdapt4->push_back({{1, 1}}); // (1,1) position in square
m_meshes.push_back(make_pair(newMeshAndPoints.first[3], iposToAdapt4));
// store new points
m_points.resize(nbpt + newMeshAndPoints.second.size());
for (size_t i = nbpt; i < nbpt + newMeshAndPoints.second.size(); ++i)
m_points[i] = newMeshAndPoints.second[i - nbpt];
}
list< pair< shared_ptr< Mesh2D>, shared_ptr<vector<ArrayXi > > > > GridAdapt2D::intersect(const double &p_xMin, const double &p_xMax, const double &p_yMin, const double &p_yMax) const
{
list< pair< shared_ptr< Mesh2D>, shared_ptr<vector<ArrayXi > > > > meshRet;
for (const auto &meshAPos : m_meshes)
{
shared_ptr< Mesh2D> mesh = meshAPos.first ;
bool bX = ((mesh->getXL() - p_xMin) * (mesh->getXL() - p_xMax) <= 0.) || ((mesh->getXR() - p_xMin) * (mesh->getXR() - p_xMax) <= 0.) || ((mesh->getXL() <= p_xMin + tiny) && (mesh->getXR() >= p_xMax - tiny));
bool bY = ((mesh->getYB() - p_yMin) * (mesh->getYB() - p_yMax) <= 0.) || ((mesh->getYT() - p_yMin) * (mesh->getYT() - p_yMax) <= 0.) || ((mesh->getYB() <= p_yMin + tiny) && (mesh->getYT() >= p_yMax - tiny));
if (bX && bY)
{
meshRet.push_back(make_pair(mesh, meshAPos.second));
}
}
return meshRet;
}
shared_ptr< Mesh2D> GridAdapt2D::getMeshWithPoint(const double &p_x, const double &p_y) const
{
for (const auto &meshAPos : m_meshes)
{
shared_ptr< Mesh2D> mesh = meshAPos.first ;
if ((mesh->getXL() <= p_x) && (mesh->getXR() >= p_x) && (mesh->getYB() <= p_y) && (mesh->getYT() >= p_y))
{
return mesh;
}
}
abort();
}
|
