Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
gmsh
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Terms and privacy
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Larry Price
gmsh
Commits
f9107096
Commit
f9107096
authored
13 years ago
by
Bastien Gorissen
Browse files
Options
Downloads
Patches
Plain Diff
Remove old smoothing code
parent
8e594ae0
No related branches found
No related tags found
No related merge requests found
Changes
3
Hide whitespace changes
Inline
Side-by-side
Showing
3 changed files
Mesh/CMakeLists.txt
+0
-1
0 additions, 1 deletion
Mesh/CMakeLists.txt
Mesh/HighOrder.cpp
+129
-105
129 additions, 105 deletions
Mesh/HighOrder.cpp
Mesh/HighOrder.h
+1
-5
1 addition, 5 deletions
Mesh/HighOrder.h
with
130 additions
and
111 deletions
Mesh/CMakeLists.txt
+
0
−
1
View file @
f9107096
...
...
@@ -24,7 +24,6 @@ CenterlineField.cpp
BoundaryLayers.cpp
BDS.cpp
HighOrder.cpp
highOrderSmoother.cpp
highOrderTools.cpp
meshPartition.cpp
meshRefine.cpp
...
...
This diff is collapsed.
Click to expand it.
Mesh/HighOrder.cpp
+
129
−
105
View file @
f9107096
...
...
@@ -8,7 +8,6 @@
//
#include
"HighOrder.h"
#include
"highOrderSmoother.h"
#include
"highOrderTools.h"
#include
"MLine.h"
#include
"MTriangle.h"
...
...
@@ -55,7 +54,7 @@ static bool computeEquidistantParameters1(GEdge *ge, double u0, double uN, int N
p0
.
y
()
+
i
*
du
*
(
p1
.
y
()
-
p0
.
y
()),
p0
.
z
()
+
i
*
du
*
(
p1
.
z
()
-
p0
.
z
()));
double
t
;
GPoint
gp
=
ge
->
closestPoint
(
pi
,
t
);
GPoint
gp
=
ge
->
closestPoint
(
pi
,
t
);
u
[
i
]
=
gp
.
u
();
// printf("going to %g %g u %g\n",pi.x(),pi.y(),gp.u());
}
...
...
@@ -257,8 +256,7 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
static
void
getEdgeVertices
(
GEdge
*
ge
,
MElement
*
ele
,
std
::
vector
<
MVertex
*>
&
ve
,
edgeContainer
&
edgeVertices
,
bool
linear
,
int
nPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
)
int
nPts
=
1
)
{
if
(
ge
->
geomType
()
==
GEntity
::
DiscreteCurve
||
ge
->
geomType
()
==
GEntity
::
BoundaryLayerCurve
)
...
...
@@ -280,7 +278,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
bool
reparamOK
=
true
;
if
(
!
linear
)
{
reparamOK
&=
reparamMeshVertexOnEdge
(
v0
,
ge
,
u0
);
if
(
ge
->
periodic
(
0
)
&&
v1
==
ge
->
getEndVertex
()
->
mesh_vertices
[
0
])
if
(
ge
->
periodic
(
0
)
&&
ge
->
getEndVertex
()
->
getNumMeshVertices
()
>
0
&&
v1
==
ge
->
getEndVertex
()
->
mesh_vertices
[
0
])
u1
=
ge
->
parBounds
(
0
).
high
();
else
reparamOK
&=
reparamMeshVertexOnEdge
(
v1
,
ge
,
u1
);
...
...
@@ -321,11 +319,6 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
GPoint
pc
=
ge
->
point
(
US
[
count
]);
v
=
new
MEdgeVertex
(
pc
.
x
(),
pc
.
y
(),
pc
.
z
(),
ge
,
pc
.
u
());
// printf("Edge %d(%g) %d(%g) new vertex %g %g at %g\n",v0->getNum(),u0,v1->getNum(),u1,v->x(),v->y(), US[count]);
if
(
displ2D
||
displ3D
){
SPoint3
pc2
=
edge
.
interpolate
(
t
);
if
(
displ2D
)
displ2D
->
add
(
v
,
SVector3
(
pc2
.
x
(),
pc2
.
y
(),
pc2
.
z
()));
if
(
displ3D
)
displ3D
->
add
(
v
,
SVector3
(
pc2
.
x
(),
pc2
.
y
(),
pc2
.
z
()));
}
}
temp
.
push_back
(
v
);
// this destroys the ordering of the mesh vertices on the edge
...
...
@@ -343,9 +336,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
static
void
getEdgeVertices
(
GFace
*
gf
,
MElement
*
ele
,
std
::
vector
<
MVertex
*>
&
ve
,
edgeContainer
&
edgeVertices
,
bool
linear
,
int
nPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
)
int
nPts
=
1
)
{
if
(
gf
->
geomType
()
==
GEntity
::
DiscreteSurface
||
gf
->
geomType
()
==
GEntity
::
BoundaryLayerSurface
)
...
...
@@ -383,11 +374,6 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
else
{
GPoint
pc
=
gf
->
point
(
US
[
j
+
1
],
VS
[
j
+
1
]);
v
=
new
MFaceVertex
(
pc
.
x
(),
pc
.
y
(),
pc
.
z
(),
gf
,
US
[
j
+
1
],
VS
[
j
+
1
]);
if
(
displ2D
||
displ3D
){
SPoint3
pc2
=
edge
.
interpolate
(
t
);
if
(
displ3D
)
displ3D
->
add
(
v
,
SVector3
(
pc2
.
x
(),
pc2
.
y
(),
pc2
.
z
()));
if
(
displ2D
)
displ2D
->
add
(
v
,
SVector3
(
pc2
.
x
(),
pc2
.
y
(),
pc2
.
z
()));
}
}
temp
.
push_back
(
v
);
gf
->
mesh_vertices
.
push_back
(
v
);
...
...
@@ -403,8 +389,7 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
static
void
getEdgeVertices
(
GRegion
*
gr
,
MElement
*
ele
,
std
::
vector
<
MVertex
*>
&
ve
,
edgeContainer
&
edgeVertices
,
bool
linear
,
int
nPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
)
int
nPts
=
1
)
{
for
(
int
i
=
0
;
i
<
ele
->
getNumEdges
();
i
++
){
MEdge
edge
=
ele
->
getEdge
(
i
);
...
...
@@ -435,8 +420,7 @@ static void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
static
void
getFaceVertices
(
GFace
*
gf
,
MElement
*
incomplete
,
MElement
*
ele
,
std
::
vector
<
MVertex
*>
&
vf
,
faceContainer
&
faceVertices
,
bool
linear
,
int
nPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
)
bool
linear
,
int
nPts
=
1
)
{
if
(
gf
->
geomType
()
==
GEntity
::
DiscreteSurface
||
gf
->
geomType
()
==
GEntity
::
BoundaryLayerSurface
)
...
...
@@ -492,11 +476,6 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
else
{
v
=
new
MVertex
(
X
,
Y
,
Z
,
gf
);
}
if
(
displ3D
||
displ2D
){
SPoint3
pc2
=
face
.
interpolate
(
t1
,
t2
);
if
(
displ3D
)
displ3D
->
add
(
v
,
SVector3
(
pc2
.
x
(),
pc2
.
y
(),
pc2
.
z
()));
if
(
displ2D
)
displ2D
->
add
(
v
,
SVector3
(
pc2
.
x
(),
pc2
.
y
(),
pc2
.
z
()));
}
}
// should be expensive -> induces a new search each time
vtcs
.
push_back
(
v
);
...
...
@@ -775,6 +754,25 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
}
start
=
(
nPts
+
2
)
*
(
nPts
+
2
)
*
(
nPts
+
2
)
-
(
nPts
)
*
(
nPts
)
*
(
nPts
)
;
break
;
case
TYPE_PYR
:
printf
(
"aaaaaand...
\n
"
);
switch
(
nPts
){
case
0
:
case
1
:
return
;
case
2
:
points
=
polynomialBases
::
find
(
MSH_PYR_30
)
->
points
;
break
;
case
3
:
points
=
polynomialBases
::
find
(
MSH_PYR_55
)
->
points
;
break
;
case
4
:
points
=
polynomialBases
::
find
(
MSH_PYR_91
)
->
points
;
break
;
case
5
:
points
=
polynomialBases
::
find
(
MSH_PYR_140
)
->
points
;
break
;
case
6
:
points
=
polynomialBases
::
find
(
MSH_PYR_204
)
->
points
;
break
;
case
7
:
points
=
polynomialBases
::
find
(
MSH_PYR_285
)
->
points
;
break
;
case
8
:
points
=
polynomialBases
::
find
(
MSH_PYR_385
)
->
points
;
break
;
default
:
Msg
::
Error
(
"getRegionVertices not implemented for order %i"
,
nPts
+
1
);
break
;
}
start
=
(
nPts
+
2
)
*
(
(
nPts
+
2
)
+
1
)
*
(
2
*
(
nPts
+
2
)
+
1
)
/
6
-
(
nPts
-
1
)
*
(
(
nPts
-
1
)
+
1
)
*
(
2
*
(
nPts
-
1
)
+
1
)
/
6
;
break
;
}
for
(
int
k
=
start
;
k
<
points
.
size1
();
k
++
){
...
...
@@ -791,14 +789,13 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
}
static
void
setHighOrder
(
GEdge
*
ge
,
edgeContainer
&
edgeVertices
,
bool
linear
,
int
nbPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
)
int
nbPts
=
1
)
{
std
::
vector
<
MLine
*>
lines2
;
for
(
unsigned
int
i
=
0
;
i
<
ge
->
lines
.
size
();
i
++
){
MLine
*
l
=
ge
->
lines
[
i
];
std
::
vector
<
MVertex
*>
ve
;
getEdgeVertices
(
ge
,
l
,
ve
,
edgeVertices
,
linear
,
nbPts
,
displ2D
,
displ3D
);
getEdgeVertices
(
ge
,
l
,
ve
,
edgeVertices
,
linear
,
nbPts
);
if
(
nbPts
==
1
)
lines2
.
push_back
(
new
MLine3
(
l
->
getVertex
(
0
),
l
->
getVertex
(
1
),
ve
[
0
]));
else
...
...
@@ -812,12 +809,10 @@ static void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
MTriangle
*
setHighOrder
(
MTriangle
*
t
,
GFace
*
gf
,
edgeContainer
&
edgeVertices
,
faceContainer
&
faceVertices
,
bool
linear
,
bool
incomplete
,
int
nPts
,
highOrderSmoother
*
displ2D
,
highOrderSmoother
*
displ3D
)
bool
linear
,
bool
incomplete
,
int
nPts
)
{
std
::
vector
<
MVertex
*>
ve
,
vf
;
getEdgeVertices
(
gf
,
t
,
ve
,
edgeVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
getEdgeVertices
(
gf
,
t
,
ve
,
edgeVertices
,
linear
,
nPts
);
if
(
nPts
==
1
){
return
new
MTriangle6
(
t
->
getVertex
(
0
),
t
->
getVertex
(
1
),
t
->
getVertex
(
2
),
ve
[
0
],
ve
[
1
],
ve
[
2
]);
...
...
@@ -828,17 +823,9 @@ MTriangle *setHighOrder(MTriangle *t, GFace *gf,
ve
,
nPts
+
1
);
}
else
{
if
(
displ2D
&&
gf
->
geomType
()
==
GEntity
::
Plane
){
MTriangle
incpl
(
t
->
getVertex
(
0
),
t
->
getVertex
(
1
),
t
->
getVertex
(
2
));
getFaceVertices
(
gf
,
&
incpl
,
t
,
vf
,
faceVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
}
else
{
MTriangleN
incpl
(
t
->
getVertex
(
0
),
t
->
getVertex
(
1
),
t
->
getVertex
(
2
),
ve
,
nPts
+
1
);
getFaceVertices
(
gf
,
&
incpl
,
t
,
vf
,
faceVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
}
MTriangleN
incpl
(
t
->
getVertex
(
0
),
t
->
getVertex
(
1
),
t
->
getVertex
(
2
),
ve
,
nPts
+
1
);
getFaceVertices
(
gf
,
&
incpl
,
t
,
vf
,
faceVertices
,
linear
,
nPts
);
ve
.
insert
(
ve
.
end
(),
vf
.
begin
(),
vf
.
end
());
return
new
MTriangleN
(
t
->
getVertex
(
0
),
t
->
getVertex
(
1
),
t
->
getVertex
(
2
),
ve
,
nPts
+
1
);
...
...
@@ -849,12 +836,10 @@ MTriangle *setHighOrder(MTriangle *t, GFace *gf,
static
MQuadrangle
*
setHighOrder
(
MQuadrangle
*
q
,
GFace
*
gf
,
edgeContainer
&
edgeVertices
,
faceContainer
&
faceVertices
,
bool
linear
,
bool
incomplete
,
int
nPts
,
highOrderSmoother
*
displ2D
,
highOrderSmoother
*
displ3D
)
bool
linear
,
bool
incomplete
,
int
nPts
)
{
std
::
vector
<
MVertex
*>
ve
,
vf
;
getEdgeVertices
(
gf
,
q
,
ve
,
edgeVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
getEdgeVertices
(
gf
,
q
,
ve
,
edgeVertices
,
linear
,
nPts
);
if
(
incomplete
){
if
(
nPts
==
1
){
return
new
MQuadrangle8
(
q
->
getVertex
(
0
),
q
->
getVertex
(
1
),
q
->
getVertex
(
2
),
...
...
@@ -866,15 +851,9 @@ static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
}
}
else
{
if
(
displ2D
&&
gf
->
geomType
()
==
GEntity
::
Plane
){
MQuadrangle
incpl
(
q
->
getVertex
(
0
),
q
->
getVertex
(
1
),
q
->
getVertex
(
2
),
q
->
getVertex
(
3
));
getFaceVertices
(
gf
,
&
incpl
,
q
,
vf
,
faceVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
}
else
{
MQuadrangleN
incpl
(
q
->
getVertex
(
0
),
q
->
getVertex
(
1
),
q
->
getVertex
(
2
),
q
->
getVertex
(
3
),
ve
,
nPts
+
1
);
getFaceVertices
(
gf
,
&
incpl
,
q
,
vf
,
faceVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
}
MQuadrangleN
incpl
(
q
->
getVertex
(
0
),
q
->
getVertex
(
1
),
q
->
getVertex
(
2
),
q
->
getVertex
(
3
),
ve
,
nPts
+
1
);
getFaceVertices
(
gf
,
&
incpl
,
q
,
vf
,
faceVertices
,
linear
,
nPts
);
ve
.
insert
(
ve
.
end
(),
vf
.
begin
(),
vf
.
end
());
if
(
nPts
==
1
){
return
new
MQuadrangle9
(
q
->
getVertex
(
0
),
q
->
getVertex
(
1
),
q
->
getVertex
(
2
),
...
...
@@ -889,14 +868,13 @@ static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
static
void
setHighOrder
(
GFace
*
gf
,
edgeContainer
&
edgeVertices
,
faceContainer
&
faceVertices
,
bool
linear
,
bool
incomplete
,
int
nPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
)
int
nPts
=
1
)
{
std
::
vector
<
MTriangle
*>
triangles2
;
for
(
unsigned
int
i
=
0
;
i
<
gf
->
triangles
.
size
();
i
++
){
MTriangle
*
t
=
gf
->
triangles
[
i
];
MTriangle
*
tNew
=
setHighOrder
(
t
,
gf
,
edgeVertices
,
faceVertices
,
linear
,
incomplete
,
nPts
,
displ2D
,
displ3D
);
nPts
);
triangles2
.
push_back
(
tNew
);
delete
t
;
}
...
...
@@ -905,7 +883,7 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
for
(
unsigned
int
i
=
0
;
i
<
gf
->
quadrangles
.
size
();
i
++
){
MQuadrangle
*
q
=
gf
->
quadrangles
[
i
];
MQuadrangle
*
qNew
=
setHighOrder
(
q
,
gf
,
edgeVertices
,
faceVertices
,
linear
,
incomplete
,
nPts
,
displ2D
,
displ3D
);
incomplete
,
nPts
);
quadrangles2
.
push_back
(
qNew
);
delete
q
;
}
...
...
@@ -915,14 +893,13 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
static
void
setHighOrder
(
GRegion
*
gr
,
edgeContainer
&
edgeVertices
,
faceContainer
&
faceVertices
,
bool
linear
,
bool
incomplete
,
int
nPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
)
int
nPts
=
1
)
{
std
::
vector
<
MTetrahedron
*>
tetrahedra2
;
for
(
unsigned
int
i
=
0
;
i
<
gr
->
tetrahedra
.
size
();
i
++
){
MTetrahedron
*
t
=
gr
->
tetrahedra
[
i
];
std
::
vector
<
MVertex
*>
ve
,
vf
,
vr
;
getEdgeVertices
(
gr
,
t
,
ve
,
edgeVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
getEdgeVertices
(
gr
,
t
,
ve
,
edgeVertices
,
linear
,
nPts
);
if
(
nPts
==
1
){
tetrahedra2
.
push_back
(
new
MTetrahedron10
(
t
->
getVertex
(
0
),
t
->
getVertex
(
1
),
t
->
getVertex
(
2
),
...
...
@@ -947,7 +924,7 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
for
(
unsigned
int
i
=
0
;
i
<
gr
->
hexahedra
.
size
();
i
++
){
MHexahedron
*
h
=
gr
->
hexahedra
[
i
];
std
::
vector
<
MVertex
*>
ve
,
vf
,
vr
;
getEdgeVertices
(
gr
,
h
,
ve
,
edgeVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
getEdgeVertices
(
gr
,
h
,
ve
,
edgeVertices
,
linear
,
nPts
);
if
(
nPts
==
1
){
if
(
incomplete
){
hexahedra2
.
push_back
...
...
@@ -990,8 +967,8 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
std
::
vector
<
MPrism
*>
prisms2
;
for
(
unsigned
int
i
=
0
;
i
<
gr
->
prisms
.
size
();
i
++
){
MPrism
*
p
=
gr
->
prisms
[
i
];
std
::
vector
<
MVertex
*>
ve
,
vf
;
getEdgeVertices
(
gr
,
p
,
ve
,
edgeVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
std
::
vector
<
MVertex
*>
ve
,
vf
,
vr
;
getEdgeVertices
(
gr
,
p
,
ve
,
edgeVertices
,
linear
,
nPts
);
if
(
incomplete
){
prisms2
.
push_back
(
new
MPrism15
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
...
...
@@ -999,12 +976,31 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
ve
[
0
],
ve
[
1
],
ve
[
2
],
ve
[
3
],
ve
[
4
],
ve
[
5
],
ve
[
6
],
ve
[
7
],
ve
[
8
]));
}
else
{
getFaceVertices
(
gr
,
p
,
vf
,
faceVertices
,
edgeVertices
,
linear
,
nPts
);
prisms2
.
push_back
(
new
MPrism18
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
p
->
getVertex
(
5
),
ve
[
0
],
ve
[
1
],
ve
[
2
],
ve
[
3
],
ve
[
4
],
ve
[
5
],
ve
[
6
],
ve
[
7
],
ve
[
8
],
vf
[
0
],
vf
[
1
],
vf
[
2
]));
if
(
nPts
==
1
)
{
getFaceVertices
(
gr
,
p
,
vf
,
faceVertices
,
edgeVertices
,
linear
,
nPts
);
prisms2
.
push_back
(
new
MPrism18
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
p
->
getVertex
(
5
),
ve
[
0
],
ve
[
1
],
ve
[
2
],
ve
[
3
],
ve
[
4
],
ve
[
5
],
ve
[
6
],
ve
[
7
],
ve
[
8
],
vf
[
0
],
vf
[
1
],
vf
[
2
]));
}
else
{
//getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
//ve.insert(ve.end(), vf.begin(), vf.end());
//MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2), p->getVertex(3),
// p->getVertex(4), p->getVertex(5), ve, nPts + 1);
//getRegionVertices(gr, &incpl, p, vr, linear, nPts);
//if (nPts == 0) {
// printf("Screwed\n");
// }
//ve.insert(ve.end(), vr.begin(), vr.end());
//MPrism* n = new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
// p->getVertex(3), p->getVertex(4), p->getVertex(5),
// ve, nPts+1);
//if (!mappingIsInvertible(n))
// Msg::Warning("Found invalid curved volume element (# %d in list)", i);
//prisms2.push_back(n);
}
}
delete
p
;
}
...
...
@@ -1013,20 +1009,35 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
std
::
vector
<
MPyramid
*>
pyramids2
;
for
(
unsigned
int
i
=
0
;
i
<
gr
->
pyramids
.
size
();
i
++
){
MPyramid
*
p
=
gr
->
pyramids
[
i
];
std
::
vector
<
MVertex
*>
ve
,
vf
;
getEdgeVertices
(
gr
,
p
,
ve
,
edgeVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
if
(
incomplete
){
pyramids2
.
push_back
(
new
MPyramid13
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
ve
[
0
],
ve
[
1
],
ve
[
2
],
ve
[
3
],
ve
[
4
],
ve
[
5
],
ve
[
6
],
ve
[
7
]));
std
::
vector
<
MVertex
*>
ve
,
vf
,
vr
;
getEdgeVertices
(
gr
,
p
,
ve
,
edgeVertices
,
linear
,
nPts
);
if
(
nPts
==
1
)
{
if
(
incomplete
){
pyramids2
.
push_back
(
new
MPyramid13
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
ve
[
0
],
ve
[
1
],
ve
[
2
],
ve
[
3
],
ve
[
4
],
ve
[
5
],
ve
[
6
],
ve
[
7
]));
}
else
{
getFaceVertices
(
gr
,
p
,
vf
,
faceVertices
,
edgeVertices
,
linear
,
nPts
);
pyramids2
.
push_back
(
new
MPyramid14
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
ve
[
0
],
ve
[
1
],
ve
[
2
],
ve
[
3
],
ve
[
4
],
ve
[
5
],
ve
[
6
],
ve
[
7
],
vf
[
0
]));
}
}
else
{
else
{
getFaceVertices
(
gr
,
p
,
vf
,
faceVertices
,
edgeVertices
,
linear
,
nPts
);
pyramids2
.
push_back
(
new
MPyramid14
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
ve
[
0
],
ve
[
1
],
ve
[
2
],
ve
[
3
],
ve
[
4
],
ve
[
5
],
ve
[
6
],
ve
[
7
],
vf
[
0
]));
ve
.
insert
(
ve
.
end
(),
vf
.
begin
(),
vf
.
end
());
MPyramidN
incpl
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
ve
,
nPts
+
1
);
getRegionVertices
(
gr
,
&
incpl
,
p
,
vr
,
linear
,
nPts
);
ve
.
insert
(
ve
.
end
(),
vr
.
begin
(),
vr
.
end
());
printf
(
"Hmmm, let's see : %d
\n
"
,
nPts
);
MPyramid
*
n
=
new
MPyramidN
(
p
->
getVertex
(
0
),
p
->
getVertex
(
1
),
p
->
getVertex
(
2
),
p
->
getVertex
(
3
),
p
->
getVertex
(
4
),
ve
,
nPts
+
1
);
pyramids2
.
push_back
(
n
);
}
delete
p
;
}
...
...
@@ -1247,25 +1258,30 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// m->writeMSH("BEFORE.msh");
highOrderSmoother
*
displ2D
=
0
;
highOrderSmoother
*
displ3D
=
0
;
/*
if(CTX::instance()->mesh.smoothInternalEdges){
displ2D = new highOrderSmoother(2);
displ3D = new highOrderSmoother(3);
}
*/
// then create new second order vertices/elements
edgeContainer
edgeVertices
;
faceContainer
faceVertices
;
Msg
::
StatusBar
(
2
,
true
,
"Meshing curves order %d..."
,
order
);
for
(
GModel
::
eiter
it
=
m
->
firstEdge
();
it
!=
m
->
lastEdge
();
++
it
)
setHighOrder
(
*
it
,
edgeVertices
,
linear
,
nPts
,
displ2D
,
displ3D
);
int
counter
=
1
;
for
(
GModel
::
eiter
it
=
m
->
firstEdge
();
it
!=
m
->
lastEdge
();
++
it
)
{
Msg
::
StatusBar
(
2
,
true
,
"Meshing curves order %d (%i/%i)..."
,
order
,
counter
,
m
->
getNumEdges
());
counter
++
;
setHighOrder
(
*
it
,
edgeVertices
,
linear
,
nPts
);
}
Msg
::
StatusBar
(
2
,
true
,
"Meshing surfaces order %d..."
,
order
);
for
(
GModel
::
fiter
it
=
m
->
firstFace
();
it
!=
m
->
lastFace
();
++
it
)
setHighOrder
(
*
it
,
edgeVertices
,
faceVertices
,
linear
,
incomplete
,
nPts
,
displ2D
,
displ3D
);
counter
=
1
;
for
(
GModel
::
fiter
it
=
m
->
firstFace
();
it
!=
m
->
lastFace
();
++
it
)
{
Msg
::
StatusBar
(
2
,
true
,
"Meshing surfaces order %d (%i/%i)..."
,
order
,
counter
,
m
->
getNumFaces
());
counter
++
;
setHighOrder
(
*
it
,
edgeVertices
,
faceVertices
,
linear
,
incomplete
,
nPts
);
}
highOrderTools
hot
(
m
);
...
...
@@ -1277,18 +1293,22 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// printJacobians(m, "smoothness_b.pos");
// m->writeMSH("RAW.msh");
if
(
displ2D
){
/*
if (displ2D){
checkHighOrderTriangles("Before optimization", m, bad, worst);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
displ2D->optimize(*it,edgeVertices,faceVertices);
checkHighOrderTriangles("After optimization", m, bad, worst);
}
}
*/
Msg
::
StatusBar
(
2
,
true
,
"Meshing volumes order %d..."
,
order
);
for
(
GModel
::
riter
it
=
m
->
firstRegion
();
it
!=
m
->
lastRegion
();
++
it
)
setHighOrder
(
*
it
,
edgeVertices
,
faceVertices
,
linear
,
incomplete
,
nPts
,
displ2D
,
displ3D
);
counter
=
1
;
for
(
GModel
::
riter
it
=
m
->
firstRegion
();
it
!=
m
->
lastRegion
();
++
it
)
{
Msg
::
StatusBar
(
2
,
true
,
"Meshing volumes order %d (%i/%i)..."
,
order
,
counter
,
m
->
getNumRegions
());
counter
++
;
setHighOrder
(
*
it
,
edgeVertices
,
faceVertices
,
linear
,
incomplete
,
nPts
);
}
/*
// smooth the 3D regions
if (displ3D){
checkHighOrderTetrahedron("Before optimization", m, bad, worst);
...
...
@@ -1296,9 +1316,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
displ3D->smooth(*it);
checkHighOrderTetrahedron("After optimization", m, bad, worst);
}
if
(
displ2D
)
delete
displ2D
;
if
(
displ3D
)
delete
displ3D
;
*/
double
t2
=
Cpu
();
...
...
@@ -1306,7 +1324,13 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// m->writeMSH("SMOOTHED.msh");
// FIXME !!
if
(
!
linear
){
//hot.ensureMinimumDistorsion(0.1);
for
(
GModel
::
fiter
it
=
m
->
firstFace
();
it
!=
m
->
lastFace
();
++
it
)
{
std
::
vector
<
MElement
*>
v
;
v
.
insert
(
v
.
begin
(),
(
*
it
)
->
triangles
.
begin
(),
(
*
it
)
->
triangles
.
end
());
v
.
insert
(
v
.
end
(),
(
*
it
)
->
quadrangles
.
begin
(),
(
*
it
)
->
quadrangles
.
end
());
hot
.
applySmoothingTo
(
v
,
-
10000000.
,
true
);
}
hot
.
ensureMinimumDistorsion
(
0.1
);
checkHighOrderTriangles
(
"Final mesh"
,
m
,
bad
,
worst
);
checkHighOrderTetrahedron
(
"Final mesh"
,
m
,
bad
,
worst
);
}
...
...
This diff is collapsed.
Click to expand it.
Mesh/HighOrder.h
+
1
−
5
View file @
f9107096
...
...
@@ -19,16 +19,12 @@ typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MVertex*> > edgeCont
// are the high order representation of the face
typedef
std
::
map
<
MFace
,
std
::
vector
<
MVertex
*>
,
Less_Face
>
faceContainer
;
class
highOrderSmoother
;
void
SetOrder1
(
GModel
*
m
);
void
SetOrderN
(
GModel
*
m
,
int
order
,
bool
linear
=
true
,
bool
incomplete
=
false
);
MTriangle
*
setHighOrder
(
MTriangle
*
t
,
GFace
*
gf
,
edgeContainer
&
edgeVertices
,
faceContainer
&
faceVertices
,
bool
linear
,
bool
incomplete
,
int
nPts
=
1
,
highOrderSmoother
*
displ2D
=
0
,
highOrderSmoother
*
displ3D
=
0
);
bool
linear
,
bool
incomplete
,
int
nPts
=
1
);
void
checkHighOrderTriangles
(
const
char
*
cc
,
GModel
*
m
,
std
::
vector
<
MElement
*>
&
bad
,
double
&
minJGlob
);
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment