Skip to content

Commit c4e9d5d

Browse files
SteveMaas1978SteveMaas1978
committed
Cylinder and tube primitives now use new multiblock mesher.
1 parent 972b063 commit c4e9d5d

File tree

5 files changed

+212
-145
lines changed

5 files changed

+212
-145
lines changed

‎MeshTools/FECylinder.cpp

+57-138
Original file line numberDiff line numberDiff line change
@@ -35,93 +35,6 @@ SOFTWARE.*/
3535

3636
extern double gain2(double x, double r, double n);
3737

38-
class FECylinderShapeModifier
39-
{
40-
public:
41-
FECylinderShapeModifier(double radius, double ratio, bool br, double gr, int ns)
42-
{
43-
m_R = radius;
44-
m_r = ratio;
45-
m_br = br;
46-
m_gr = gr;
47-
m_ns = ns;
48-
}
49-
50-
vec3d Apply(const vec3d& r)
51-
{
52-
vec3d rn = r;
53-
54-
double R1 = m_R;
55-
double R0 = m_r;
56-
if (R0 < 0) R0 = 0;
57-
if (R0 > 1) R0 = 1;
58-
R0 *= R1;
59-
60-
double d0 = R0 / sqrt(2.0);
61-
double d1 = R1 / sqrt(2.0);
62-
63-
// project the nodes onto a cylinder
64-
vec3d r0, r1;
65-
66-
// get the nodal coordinate in the template
67-
double x = rn.x;
68-
double y = rn.y;
69-
70-
// get the max-distance
71-
double D = fmax(fabs(x), fabs(y));
72-
73-
if (D <= 1)
74-
{
75-
rn.x *= d0;
76-
rn.y *= d0;
77-
}
78-
else
79-
{
80-
// "normalize" the coordinates
81-
// with respect to the max distance
82-
double r = x / D;
83-
double s = y / D;
84-
85-
vec3d r0;
86-
if (fabs(x) >= fabs(y))
87-
{
88-
double u = x / fabs(x);
89-
r0.x = u * R1*cos(PI*0.25*s);
90-
r0.y = R1 * sin(PI*0.25*s);
91-
}
92-
else
93-
{
94-
double u = y / fabs(y);
95-
r0.y = u * R1*cos(PI*0.25*r);
96-
r0.x = R1 * sin(PI*0.25*r);
97-
}
98-
99-
vec3d r1(r*d0, s*d0, 0);
100-
double a = D - 1;
101-
102-
if (m_br)
103-
{
104-
if (a <= 0.5)
105-
a = 0.5*gain2(2 * a, m_gr, m_ns);
106-
else
107-
a = 1 - 0.5*gain2(2 - 2 * a, m_gr, m_ns);
108-
}
109-
else a = gain2(a, m_gr, m_ns);
110-
111-
rn.x = r0.x*a + r1.x*(1 - a);
112-
rn.y = r0.y*a + r1.y*(1 - a);
113-
}
114-
115-
return rn;
116-
}
117-
118-
private:
119-
double m_R, m_r;
120-
bool m_br;
121-
double m_gr;
122-
int m_ns;
123-
};
124-
12538
//-----------------------------------------------------------------------------
12639
// Constructor
12740
FECylinder::FECylinder(GCylinder* po)
@@ -192,7 +105,7 @@ FEMesh* FECylinder::BuildButterfly()
192105
int ns = m_ns;
193106
int nz = m_nz;
194107
double fz = m_gz;
195-
double fr = 1;
108+
double fr = m_gr;
196109

197110
// check parameters
198111
if (nd < 1) nd = 1;
@@ -203,45 +116,49 @@ FEMesh* FECylinder::BuildButterfly()
203116
double R1 = param.GetFloatValue(GCylinder::RADIUS);
204117
m_r = GetFloatValue(RATIO);
205118

119+
double b = R1 * sqrt(2.0) / 2.0;
120+
double a = m_r * b;
121+
double c = R1;
122+
206123
// create the MB nodes
207124
m_MBNode.resize(34);
208-
m_MBNode[ 0].m_r = vec3d( -1, -1, 0);
209-
m_MBNode[ 1].m_r = vec3d( 0, -1, 0);
210-
m_MBNode[ 2].m_r = vec3d( 1, -1, 0);
211-
m_MBNode[ 3].m_r = vec3d( -1, 0, 0);
125+
m_MBNode[ 0].m_r = vec3d( -a, -a, 0);
126+
m_MBNode[ 1].m_r = vec3d( 0, -a, 0);
127+
m_MBNode[ 2].m_r = vec3d( a, -a, 0);
128+
m_MBNode[ 3].m_r = vec3d( -a, 0, 0);
212129
m_MBNode[ 4].m_r = vec3d( 0, 0, 0);
213-
m_MBNode[ 5].m_r = vec3d( 1, 0, 0);
214-
m_MBNode[ 6].m_r = vec3d( -1, 1, 0);
215-
m_MBNode[ 7].m_r = vec3d( 0, 1, 0);
216-
m_MBNode[ 8].m_r = vec3d( 1, 1, 0);
217-
218-
m_MBNode[ 9].m_r = vec3d(-1, -1, h);
219-
m_MBNode[10].m_r = vec3d( 0, -1, h);
220-
m_MBNode[11].m_r = vec3d( 1, -1, h);
221-
m_MBNode[12].m_r = vec3d(-1, 0, h);
130+
m_MBNode[ 5].m_r = vec3d( a, 0, 0);
131+
m_MBNode[ 6].m_r = vec3d( -a, a, 0);
132+
m_MBNode[ 7].m_r = vec3d( 0, a, 0);
133+
m_MBNode[ 8].m_r = vec3d( a, a, 0);
134+
135+
m_MBNode[ 9].m_r = vec3d(-a, -a, h);
136+
m_MBNode[10].m_r = vec3d( 0, -a, h);
137+
m_MBNode[11].m_r = vec3d( a, -a, h);
138+
m_MBNode[12].m_r = vec3d(-a, 0, h);
222139
m_MBNode[13].m_r = vec3d( 0, 0, h);
223-
m_MBNode[14].m_r = vec3d( 1, 0, h);
224-
m_MBNode[15].m_r = vec3d(-1, 1, h);
225-
m_MBNode[16].m_r = vec3d( 0, 1, h);
226-
m_MBNode[17].m_r = vec3d( 1, 1, h);
227-
228-
m_MBNode[18].m_r = vec3d(-2, -2, 0);
229-
m_MBNode[19].m_r = vec3d( 0, -2, 0);
230-
m_MBNode[20].m_r = vec3d( 2, -2, 0);
231-
m_MBNode[21].m_r = vec3d( 2, 0, 0);
232-
m_MBNode[22].m_r = vec3d( 2, 2, 0);
233-
m_MBNode[23].m_r = vec3d( 0, 2, 0);
234-
m_MBNode[24].m_r = vec3d(-2, 2, 0);
235-
m_MBNode[25].m_r = vec3d(-2, 0, 0);
236-
237-
m_MBNode[26].m_r = vec3d(-2, -2, h);
238-
m_MBNode[27].m_r = vec3d( 0, -2, h);
239-
m_MBNode[28].m_r = vec3d( 2, -2, h);
240-
m_MBNode[29].m_r = vec3d( 2, 0, h);
241-
m_MBNode[30].m_r = vec3d( 2, 2, h);
242-
m_MBNode[31].m_r = vec3d( 0, 2, h);
243-
m_MBNode[32].m_r = vec3d(-2, 2, h);
244-
m_MBNode[33].m_r = vec3d(-2, 0, h);
140+
m_MBNode[14].m_r = vec3d( a, 0, h);
141+
m_MBNode[15].m_r = vec3d(-a, a, h);
142+
m_MBNode[16].m_r = vec3d( 0, a, h);
143+
m_MBNode[17].m_r = vec3d( a, a, h);
144+
145+
m_MBNode[18].m_r = vec3d(-b, -b, 0);
146+
m_MBNode[19].m_r = vec3d( 0, -c, 0);
147+
m_MBNode[20].m_r = vec3d( b, -b, 0);
148+
m_MBNode[21].m_r = vec3d( c, 0, 0);
149+
m_MBNode[22].m_r = vec3d( b, b, 0);
150+
m_MBNode[23].m_r = vec3d( 0, c, 0);
151+
m_MBNode[24].m_r = vec3d(-b, b, 0);
152+
m_MBNode[25].m_r = vec3d(-c, 0, 0);
153+
154+
m_MBNode[26].m_r = vec3d(-b, -b, h);
155+
m_MBNode[27].m_r = vec3d( 0, -c, h);
156+
m_MBNode[28].m_r = vec3d( b, -b, h);
157+
m_MBNode[29].m_r = vec3d( c, 0, h);
158+
m_MBNode[30].m_r = vec3d( b, b, h);
159+
m_MBNode[31].m_r = vec3d( 0, c, h);
160+
m_MBNode[32].m_r = vec3d(-b, b, h);
161+
m_MBNode[33].m_r = vec3d(-c, 0, h);
245162

246163
// create the MB blocks
247164
m_MBlock.resize(12);
@@ -358,6 +275,23 @@ FEMesh* FECylinder::BuildButterfly()
358275
MBFace& F7 = GetBlockFace( 5, 1); SetFaceEdgeID(F7, 3, -1, 7, 11);
359276
MBFace& F8 = GetBlockFace( 6, 1); SetFaceEdgeID(F8, 3, 8, 7, -1);
360277

278+
GetFaceEdge(F1, 0).edge.m_ntype = EDGE_ZARC;
279+
GetFaceEdge(F1, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F1, 2).m_winding = -1;
280+
GetFaceEdge(F2, 0).edge.m_ntype = EDGE_ZARC;
281+
GetFaceEdge(F2, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F2, 2).m_winding = -1;
282+
GetFaceEdge(F3, 0).edge.m_ntype = EDGE_ZARC;
283+
GetFaceEdge(F3, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F3, 2).m_winding = -1;
284+
GetFaceEdge(F4, 0).edge.m_ntype = EDGE_ZARC;
285+
GetFaceEdge(F4, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F4, 2).m_winding = -1;
286+
GetFaceEdge(F5, 0).edge.m_ntype = EDGE_ZARC;
287+
GetFaceEdge(F5, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F5, 2).m_winding = -1;
288+
GetFaceEdge(F6, 0).edge.m_ntype = EDGE_ZARC;
289+
GetFaceEdge(F6, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F6, 2).m_winding = -1;
290+
GetFaceEdge(F7, 0).edge.m_ntype = EDGE_ZARC;
291+
GetFaceEdge(F7, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F7, 2).m_winding = -1;
292+
GetFaceEdge(F8, 0).edge.m_ntype = EDGE_ZARC;
293+
GetFaceEdge(F8, 2).edge.m_ntype = EDGE_ZARC; GetFaceEdge(F8, 2).m_winding = -1;
294+
361295
m_MBNode[21].SetID(0);
362296
m_MBNode[23].SetID(1);
363297
m_MBNode[25].SetID(2);
@@ -370,27 +304,12 @@ FEMesh* FECylinder::BuildButterfly()
370304
// create the MB
371305
FEMesh* pm = FEMultiBlockMesh::BuildMesh();
372306

373-
// apply a global shape modifier
374-
FECylinderShapeModifier* mod = new FECylinderShapeModifier(R1, m_r, m_br, m_gr, m_ns);
375-
for (int i = 0; i < pm->Nodes(); ++i)
376-
{
377-
vec3d r0 = pm->Node(i).pos();
378-
vec3d rn = mod->Apply(r0);
379-
pm->Node(i).pos(rn);
380-
}
381-
382-
// update the mesh
383-
pm->UpdateMesh();
384-
385307
// the Multi-block mesher will assign a different smoothing ID
386308
// to each face, but we don't want that here.
387309
// For now, we autosmooth the mesh although we should think of a
388310
// better way
389311
pm->AutoSmooth(60);
390312

391-
// cleanup
392-
delete mod;
393-
394313
return pm;
395314
}
396315

‎MeshTools/FEMultiBlockMesh.cpp

+23-1
Original file line numberDiff line numberDiff line change
@@ -42,14 +42,15 @@ void MBBlock::SetNodes(int n1,int n2,int n3,int n4,int n5,int n6,int n7,int n8)
4242
m_node[7] = n8;
4343
}
4444

45-
void MBBlock::SetZoning(double gx, double gy, double gz, bool bx, bool by, bool bz)
45+
MBBlock& MBBlock::SetZoning(double gx, double gy, double gz, bool bx, bool by, bool bz)
4646
{
4747
m_gx = gx;
4848
m_gy = gy;
4949
m_gz = gz;
5050
m_bx = bx;
5151
m_by = by;
5252
m_bz = bz;
53+
return *this;
5354
}
5455

5556
//-----------------------------------------------------------------------------
@@ -185,6 +186,18 @@ void FEMultiBlockMesh::BuildNodes(FEMesh *pm)
185186
case EDGE_LINE:
186187
pn->r = r1 * (1 - r) + r2 * r;
187188
break;
189+
case EDGE_ZARC:
190+
{
191+
vec2d c(0, 0);
192+
vec2d a(r1.x, r1.y);
193+
vec2d b(r2.x, r2.y);
194+
195+
// create an arc object
196+
GM_CIRCLE_ARC ca(c, a, b, e.m_winding);
197+
vec2d p = ca.Point(r);
198+
pn->r = vec3d(p.x, p.y, r1.z);
199+
}
200+
break;
188201
case EDGE_3P_CIRC_ARC:
189202
{
190203
vec3d r0 = m_MBNode[e.edge.m_cnode].m_r;
@@ -1228,6 +1241,15 @@ MBNode& FEMultiBlockMesh::AddNode(const vec3d& r, int nodeType)
12281241
return m_MBNode[m_MBNode.size() - 1];
12291242
}
12301243

1244+
//-----------------------------------------------------------------------------
1245+
MBBlock& FEMultiBlockMesh::AddBlock(int n0, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
1246+
{
1247+
MBBlock b;
1248+
b.SetNodes(n0, n1, n2, n3, n4, n5, n6, n7);
1249+
m_MBlock.push_back(b);
1250+
return m_MBlock[m_MBlock.size() - 1];
1251+
}
1252+
12311253
//-----------------------------------------------------------------------------
12321254
MBEdge& FEMultiBlockMesh::GetEdge(int nedge)
12331255
{

‎MeshTools/FEMultiBlockMesh.h

+10-6
Original file line numberDiff line numberDiff line change
@@ -164,8 +164,8 @@ class MBBlock : public MBItem
164164
}
165165

166166
void SetNodes(int n1,int n2,int n3,int n4,int n5,int n6,int n7,int n8);
167-
void SetSizes(int nx, int ny, int nz) { m_nx = nx; m_ny = ny; m_nz = nz; }
168-
void SetZoning(double gx, double gy, double gz, bool bx, bool by, bool bz);
167+
MBBlock& SetSizes(int nx, int ny, int nz) { m_nx = nx; m_ny = ny; m_nz = nz; return *this; }
168+
MBBlock& SetZoning(double gx, double gy, double gz, bool bx, bool by, bool bz);
169169

170170
public:
171171
int m_node[8]; // the eight nodes of the block
@@ -192,7 +192,10 @@ class FEMultiBlockMesh : public FEMesher
192192
// build the mesh
193193
FEMesh* BuildMesh();
194194

195-
protected:
195+
MBNode& AddNode(const vec3d& r, int nodeType = NODE_VERTEX);
196+
197+
MBBlock& AddBlock(int n0, int n1, int n2, int n3, int n4, int n5, int n6, int n7);
198+
196199
// update the Multi-Block data
197200
void UpdateMB();
198201

@@ -202,7 +205,10 @@ class FEMultiBlockMesh : public FEMesher
202205

203206
MBEdge& GetEdge(int nedge);
204207

205-
MBNode& AddNode(const vec3d& r, int nodeType = NODE_VERTEX);
208+
MBBlock& GetBlock(int i) { return m_MBlock[i]; }
209+
210+
void SetBlockFaceID(MBBlock& b, int n0, int n1, int n2, int n3, int n4, int n5);
211+
void SetFaceEdgeID(MBFace& f, int n0, int n1, int n2, int n3);
206212

207213
protected:
208214
void FindBlockNeighbours();
@@ -226,8 +232,6 @@ class FEMultiBlockMesh : public FEMesher
226232
int GetBlockNodeIndex(MBBlock& b, int i, int j, int k);
227233
int GetFaceNodeIndex(MBFace& f, int i, int j);
228234
int GetEdgeNodeIndex(MBEdge& e, int i);
229-
void SetBlockFaceID(MBBlock& b, int n0, int n1, int n2, int n3, int n4, int n5);
230-
void SetFaceEdgeID(MBFace& f, int n0, int n1, int n2, int n3);
231235

232236
int GetBlockFaceNodeIndex(MBBlock& b, int nf, int i, int j);
233237
int GetFaceEdgeNodeIndex(MBFace& f, int ne, int i);

0 commit comments

Comments
 (0)