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
5e667aa6
Commit
5e667aa6
authored
14 years ago
by
Emilie Marchandise
Browse files
Options
Downloads
Patches
Plain Diff
Compounds with nonLinear conformal map
parent
cc59c00c
No related branches found
Branches containing commit
No related tags found
Tags containing commit
No related merge requests found
Changes
4
Hide whitespace changes
Inline
Side-by-side
Showing
4 changed files
Geo/GFaceCompound.cpp
+106
-75
106 additions, 75 deletions
Geo/GFaceCompound.cpp
Geo/GFaceCompound.h
+5
-0
5 additions, 0 deletions
Geo/GFaceCompound.h
Solver/eigenSolver.cpp
+1
-1
1 addition, 1 deletion
Solver/eigenSolver.cpp
benchmarks/stl/artery.geo
+1
-3
1 addition, 3 deletions
benchmarks/stl/artery.geo
with
113 additions
and
79 deletions
Geo/GFaceCompound.cpp
+
106
−
75
View file @
5e667aa6
...
@@ -660,7 +660,6 @@ bool GFaceCompound::parametrize() const
...
@@ -660,7 +660,6 @@ bool GFaceCompound::parametrize() const
// if (!noOverlap) {
// if (!noOverlap) {
// Msg::Warning("!!! Overlap: parametrization switched to 'nonLIN conformal' map");
// Msg::Warning("!!! Overlap: parametrization switched to 'nonLIN conformal' map");
// noOverlap = parametrize_conformal_nonLinear() ;
// noOverlap = parametrize_conformal_nonLinear() ;
// exit(1);
// }
// }
if
(
!
noOverlap
||
!
checkOrientation
(
0
)
){
if
(
!
noOverlap
||
!
checkOrientation
(
0
)
){
Msg
::
Warning
(
"$$$ Overlap: parametrization switched to 'harmonic' map"
);
Msg
::
Warning
(
"$$$ Overlap: parametrization switched to 'harmonic' map"
);
...
@@ -1195,18 +1194,64 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
...
@@ -1195,18 +1194,64 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
_lsys
->
clear
();
_lsys
->
clear
();
}
}
void
GFaceCompound
::
computeThetaDerivatives
(
MVertex
*
prev
,
MVertex
*
curr
,
MVertex
*
next
,
double
&
dTdu1
,
double
&
dTdv1
,
double
&
dTdu2
,
double
&
dTdv2
,
double
&
dTdu3
,
double
&
dTdv3
)
const
{
SPoint2
p1
=
getCoordinates
(
prev
);
SPoint2
p2
=
getCoordinates
(
curr
);
SPoint2
p3
=
getCoordinates
(
next
);
SVector3
va
(
p2
.
x
()
-
p1
.
x
(),
p2
.
y
()
-
p1
.
y
(),
0.
);
SVector3
vb
(
p3
.
x
()
-
p2
.
x
(),
p3
.
y
()
-
p2
.
y
(),
0.
);
SVector3
vc
(
p1
.
x
()
-
p3
.
x
(),
p1
.
y
()
-
p3
.
y
(),
0.
);
double
a
=
norm
(
va
),
b
=
norm
(
vb
),
c
=
norm
(
vc
);
SVector3
n
=
crossprod
(
va
,
vb
);
double
sign
=
1.
;
if
(
n
.
z
()
<
0.0
)
sign
=
-
1.
;
//- compute grad_uv theta
//dTheta/du1 sign*acos((a^2+b^2 -c^2)/(2*a*b))
//a=sqrt((u2-u1)^2+(v2-v1)^2))
//b=sqrt((u3-u2)^2+(v3-v2)^2))
//c=sqrt((u1-u3)^2+(v1-v3)^2))
double
u1
=
p1
.
x
();
double
v1
=
p1
.
y
();
double
u2
=
p2
.
x
();
double
v2
=
p2
.
y
();
double
u3
=
p3
.
x
();
double
v3
=
p3
.
y
();
dTdu1
=
sign
*
(
v1
*
v3
*
u1
-
v1
*
v3
*
u2
-
v3
*
v2
*
u1
+
v3
*
v2
*
u2
-
v2
*
v1
*
u1
+
2
*
u3
*
v2
*
v1
+
u2
*
SQU
(
v1
)
-
u3
*
SQU
(
v2
)
-
u3
*
SQU
(
v1
)
+
SQU
(
v2
)
*
u1
-
v2
*
v1
*
u2
)
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
sqrt
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
))
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
dTdv1
=
sign
*
(
v1
*
u1
*
u3
-
SQU
(
u1
)
*
v3
+
v1
*
SQU
(
u2
)
-
u2
*
u3
*
v1
-
v1
*
u1
*
u2
-
u2
*
u1
*
v2
+
u3
*
u2
*
v2
-
u1
*
u3
*
v2
+
2
*
v3
*
u2
*
u1
-
v3
*
SQU
(
u2
)
+
v2
*
SQU
(
u1
))
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
sqrt
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
))
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
dTdu2
=
sign
*
(
u1
*
SQU
(
v1
)
*
SQU
(
v2
)
+
u1
*
SQU
(
v1
)
*
SQU
(
v3
)
+
v2
*
CUB
(
v1
)
*
u3
-
u2
*
v2
*
CUB
(
v1
)
-
v3
*
CUB
(
v2
)
*
u3
+
CUB
(
u1
)
*
SQU
(
v3
)
+
CUB
(
u1
)
*
SQU
(
v2
)
+
5
*
u2
*
v2
*
v1
*
SQU
(
u3
)
-
3
*
SQU
(
u2
)
*
v2
*
v1
*
u3
+
u2
*
v2
*
v1
*
SQU
(
v3
)
-
2
*
u2
*
SQU
(
v2
)
*
v1
*
v3
+
u1
*
v2
*
v1
*
SQU
(
u3
)
+
3
*
u1
*
v2
*
v1
*
SQU
(
u2
)
+
u1
*
v2
*
v1
*
SQU
(
v3
)
+
u1
*
SQU
(
v2
)
*
v1
*
v3
+
4
*
u3
*
u2
*
u1
*
SQU
(
v2
)
-
u2
*
v3
*
v2
*
SQU
(
u3
)
+
3
*
SQU
(
u2
)
*
v3
*
v2
*
u3
+
u1
*
v3
*
v2
*
SQU
(
u3
)
-
3
*
u1
*
v3
*
v2
*
SQU
(
u2
)
-
u1
*
v1
*
v3
*
SQU
(
u3
)
+
v2
*
v1
*
u3
*
SQU
(
u1
)
-
u2
*
v2
*
v1
*
SQU
(
u1
)
+
v3
*
SQU
(
v2
)
*
u3
*
v1
+
v3
*
v2
*
u3
*
SQU
(
v1
)
-
3
*
u3
*
v1
*
v3
*
SQU
(
u2
)
-
u3
*
v1
*
v3
*
SQU
(
u1
)
+
u2
*
v1
*
v3
*
SQU
(
u1
)
-
4
*
u1
*
v2
*
v1
*
u3
*
u2
-
u1
*
CUB
(
v2
)
*
v1
-
u2
*
CUB
(
v3
)
*
v2
+
u2
*
SQU
(
v3
)
*
SQU
(
v2
)
+
u1
*
CUB
(
v3
)
*
v2
-
2
*
u1
*
SQU
(
v3
)
*
SQU
(
v2
)
+
u1
*
v3
*
CUB
(
v2
)
-
u1
*
v1
*
CUB
(
v3
)
-
SQU
(
v2
)
*
u3
*
SQU
(
u1
)
+
CUB
(
v2
)
*
u3
*
v1
-
2
*
SQU
(
v2
)
*
u3
*
SQU
(
v1
)
+
3
*
u3
*
SQU
(
u2
)
*
SQU
(
v1
)
+
4
*
u1
*
v1
*
v3
*
u3
*
u2
-
3
*
u1
*
v1
*
v3
*
SQU
(
u2
)
+
SQU
(
u1
)
*
u3
*
v3
*
v2
+
u2
*
SQU
(
v1
)
*
v3
*
v2
+
5
*
u2
*
SQU
(
u1
)
*
v3
*
v2
-
2
*
u1
*
SQU
(
v1
)
*
v3
*
v2
-
4
*
u3
*
u2
*
u1
*
v3
*
v2
-
2
*
u2
*
SQU
(
v2
)
*
SQU
(
u1
)
+
u2
*
SQU
(
v2
)
*
SQU
(
v1
)
-
2
*
SQU
(
u3
)
*
u2
*
SQU
(
v2
)
-
3
*
SQU
(
u3
)
*
u2
*
SQU
(
v1
)
-
u1
*
SQU
(
u3
)
*
SQU
(
v2
)
+
3
*
SQU
(
u2
)
*
u1
*
SQU
(
v3
)
+
SQU
(
v1
)
*
CUB
(
u3
)
-
2
*
u2
*
SQU
(
v1
)
*
SQU
(
v3
)
-
3
*
u2
*
SQU
(
u1
)
*
SQU
(
v3
)
+
u2
*
v1
*
v3
*
SQU
(
u3
)
-
2
*
u3
*
v2
*
v1
*
SQU
(
v3
)
-
2
*
CUB
(
u1
)
*
v3
*
v2
-
CUB
(
u2
)
*
SQU
(
v3
)
-
SQU
(
v1
)
*
CUB
(
u2
)
-
u3
*
CUB
(
v1
)
*
v3
+
2
*
CUB
(
u2
)
*
v1
*
v3
+
u2
*
CUB
(
v1
)
*
v3
+
u2
*
v1
*
CUB
(
v3
)
-
2
*
v2
*
v1
*
CUB
(
u3
)
+
u3
*
SQU
(
v1
)
*
SQU
(
v3
)
+
u3
*
SQU
(
v2
)
*
SQU
(
v3
)
+
SQU
(
v2
)
*
CUB
(
u3
))
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
dTdv2
=
-
sign
*
(
CUB
(
u3
)
*
u2
*
v2
-
u1
*
CUB
(
u3
)
*
v2
-
2
*
u1
*
u3
*
CUB
(
v2
)
+
3
*
u1
*
u3
*
SQU
(
v2
)
*
v1
-
u1
*
u3
*
v2
*
SQU
(
v1
)
+
v3
*
u1
*
u3
*
SQU
(
v1
)
-
4
*
v3
*
u1
*
u3
*
v2
*
v1
+
v1
*
u1
*
u3
*
SQU
(
v3
)
+
u3
*
u2
*
v2
*
SQU
(
v3
)
-
u2
*
u3
*
v1
*
SQU
(
v3
)
+
v3
*
CUB
(
u1
)
*
u3
-
u2
*
CUB
(
u3
)
*
v1
+
u2
*
CUB
(
u1
)
*
v2
+
v1
*
u1
*
CUB
(
u3
)
+
v1
*
SQU
(
u2
)
*
SQU
(
v3
)
+
v1
*
CUB
(
u2
)
*
u1
-
v1
*
CUB
(
u2
)
*
u3
+
2
*
v1
*
SQU
(
u2
)
*
SQU
(
u3
)
-
CUB
(
u2
)
*
u1
*
v3
-
v2
*
SQU
(
u1
)
*
u3
*
u2
+
2
*
v2
*
SQU
(
u2
)
*
u1
*
u3
-
v2
*
u2
*
u1
*
SQU
(
u3
)
-
u1
*
u3
*
v2
*
SQU
(
v3
)
+
3
*
u1
*
u3
*
SQU
(
v2
)
*
v3
+
2
*
SQU
(
u2
)
*
SQU
(
u1
)
*
v3
-
u2
*
CUB
(
u1
)
*
v3
+
4
*
u2
*
u1
*
v3
*
v2
*
v1
-
CUB
(
v1
)
*
SQU
(
u3
)
-
CUB
(
v1
)
*
SQU
(
u2
)
+
2
*
CUB
(
v1
)
*
u3
*
u2
+
2
*
v2
*
SQU
(
v1
)
*
SQU
(
u2
)
+
3
*
v2
*
SQU
(
v1
)
*
SQU
(
u3
)
-
3
*
v1
*
SQU
(
v2
)
*
SQU
(
u3
)
-
v1
*
SQU
(
u1
)
*
SQU
(
u2
)
-
v1
*
SQU
(
u1
)
*
SQU
(
u3
)
+
3
*
SQU
(
u1
)
*
SQU
(
v3
)
*
v2
+
2
*
u2
*
u1
*
CUB
(
v3
)
+
2
*
SQU
(
u2
)
*
SQU
(
v3
)
*
v2
+
v3
*
SQU
(
v1
)
*
SQU
(
u2
)
-
3
*
v3
*
SQU
(
u1
)
*
SQU
(
v2
)
-
5
*
u2
*
u1
*
SQU
(
v3
)
*
v2
+
u2
*
u1
*
v2
*
SQU
(
v1
)
-
SQU
(
u1
)
*
CUB
(
v3
)
-
SQU
(
u2
)
*
CUB
(
v3
)
+
4
*
v3
*
v2
*
v1
*
u3
*
u2
-
v3
*
SQU
(
u1
)
*
SQU
(
u3
)
+
v3
*
CUB
(
u2
)
*
u3
-
v3
*
SQU
(
u2
)
*
SQU
(
u3
)
-
SQU
(
u2
)
*
v2
*
SQU
(
u3
)
+
3
*
u2
*
u1
*
v3
*
SQU
(
v2
)
-
u2
*
u1
*
v3
*
SQU
(
v1
)
-
4
*
v3
*
v2
*
v1
*
SQU
(
u2
)
-
5
*
v2
*
SQU
(
v1
)
*
u3
*
u2
+
3
*
v1
*
SQU
(
v2
)
*
u3
*
u2
-
CUB
(
u1
)
*
u3
*
v2
-
v2
*
SQU
(
u1
)
*
SQU
(
u2
)
+
2
*
v2
*
SQU
(
u1
)
*
SQU
(
u3
)
+
SQU
(
u1
)
*
CUB
(
v2
)
+
CUB
(
v2
)
*
SQU
(
u3
)
+
2
*
v1
*
SQU
(
u1
)
*
u3
*
u2
-
3
*
v1
*
u2
*
u1
*
SQU
(
v2
)
-
v1
*
u2
*
u1
*
SQU
(
v3
)
-
v1
*
SQU
(
u2
)
*
u1
*
u3
-
v1
*
u2
*
u1
*
SQU
(
u3
)
-
v3
*
SQU
(
v1
)
*
u3
*
u2
-
3
*
v3
*
SQU
(
v2
)
*
u3
*
u2
-
v3
*
SQU
(
u1
)
*
u3
*
u2
-
v3
*
SQU
(
u2
)
*
u1
*
u3
+
2
*
v3
*
u2
*
u1
*
SQU
(
u3
))
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
dTdu3
=
-
sign
*
(
-
SQU
(
v2
)
*
u3
-
2
*
u1
*
v3
*
v2
-
u2
*
SQU
(
v3
)
+
u2
*
v3
*
v2
+
v2
*
v1
*
u3
+
u2
*
v1
*
v3
+
u1
*
SQU
(
v3
)
-
u2
*
v2
*
v1
+
u1
*
SQU
(
v2
)
+
v3
*
v2
*
u3
-
u3
*
v1
*
v3
)
/
sqrt
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
dTdv3
=
sign
*
(
v3
*
u1
*
u3
-
v1
*
SQU
(
u3
)
+
2
*
u2
*
u3
*
v1
-
v1
*
SQU
(
u2
)
+
u2
*
u1
*
v2
-
u2
*
u1
*
v3
-
u3
*
u2
*
v2
-
u3
*
u2
*
v3
-
u1
*
u3
*
v2
+
SQU
(
u2
)
*
v3
+
v2
*
SQU
(
u3
))
/
sqrt
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
}
bool
GFaceCompound
::
parametrize_conformal_nonLinear
()
const
bool
GFaceCompound
::
parametrize_conformal_nonLinear
()
const
{
{
printStuff
(
100
);
bool
converged
=
false
;
//exit(1);
//--create dofManager
//--create dofManager
dofManager
<
double
>
myAssembler
(
_lsys
);
linearSystem
<
double
>*
lsysNL
;
lsysNL
=
new
linearSystemFull
<
double
>
;
//lsysNL = new linearSystemPETSc<double>;
//lsysNL = new linearSystemCSRTaucs<double>;
dofManager
<
double
>
myAssembler
(
lsysNL
);
//--- first compute mapping harmonic
//--- first compute mapping harmonic
parametrize
(
ITERU
,
HARMONIC
);
parametrize
(
ITERU
,
HARMONIC
);
parametrize
(
ITERV
,
HARMONIC
);
parametrize
(
ITERV
,
HARMONIC
);
printStuff
(
100
);
//---order boundary vertices
//---order boundary vertices
std
::
vector
<
MVertex
*>
ordered
;
std
::
vector
<
MVertex
*>
ordered
;
...
@@ -1228,19 +1273,15 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
...
@@ -1228,19 +1273,15 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
myAssembler
.
numberVertex
(
v
,
0
,
1
);
myAssembler
.
numberVertex
(
v
,
0
,
1
);
myAssembler
.
numberVertex
(
v
,
0
,
2
);
myAssembler
.
numberVertex
(
v
,
0
,
2
);
}
}
for
(
std
::
set
<
MVertex
*>::
iterator
itv
=
fillNodes
.
begin
();
itv
!=
fillNodes
.
end
()
;
++
itv
){
MVertex
*
v
=
*
itv
;
myAssembler
.
numberVertex
(
v
,
0
,
1
);
myAssembler
.
numberVertex
(
v
,
0
,
2
);
}
MVertex
*
lag
=
new
MVertex
(
0.
,
0.
,
0.
);
MVertex
*
lag
=
new
MVertex
(
0.
,
0.
,
0.
);
myAssembler
.
numberVertex
(
lag
,
0
,
3
);
//ghost vertex for lagrange multiplier
myAssembler
.
numberVertex
(
lag
,
0
,
3
);
//ghost vertex for lagrange multiplier
//--- newton Loop
//--- newton Loop
int
nbNewton
=
3
;
int
nbNewton
=
5
;
double
lambda
=
1.e10
;
double
lambda
=
1.e5
;
double
fac
=
1.e7
;
for
(
int
iNewton
=
0
;
iNewton
<
nbNewton
;
iNewton
++
){
for
(
int
iNewton
=
0
;
iNewton
<
nbNewton
;
iNewton
++
){
//-- assemble conformal matrix
//-- assemble conformal matrix
std
::
vector
<
MElement
*>
allElems
;
std
::
vector
<
MElement
*>
allElems
;
laplaceTerm
laplace1
(
model
(),
1
,
ONE
,
&
coordinates
);
laplaceTerm
laplace1
(
model
(),
1
,
ONE
,
&
coordinates
);
...
@@ -1258,14 +1299,6 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
...
@@ -1258,14 +1299,6 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
allElems
.
push_back
((
MElement
*
)((
*
it
)
->
triangles
[
i
]));
allElems
.
push_back
((
MElement
*
)((
*
it
)
->
triangles
[
i
]));
}
}
}
}
for
(
std
::
list
<
MTriangle
*>::
iterator
it2
=
fillTris
.
begin
();
it2
!=
fillTris
.
end
();
it2
++
){
SElement
se
((
*
it2
));
laplace1
.
addToMatrix
(
myAssembler
,
&
se
);
laplace2
.
addToMatrix
(
myAssembler
,
&
se
);
cross12
.
addToMatrix
(
myAssembler
,
&
se
);
cross21
.
addToMatrix
(
myAssembler
,
&
se
);
allElems
.
push_back
((
MElement
*
)(
*
it2
));
}
groupOfElements
gr
(
allElems
);
groupOfElements
gr
(
allElems
);
laplace1
.
addToRightHandSide
(
myAssembler
,
gr
);
laplace1
.
addToRightHandSide
(
myAssembler
,
gr
);
...
@@ -1288,11 +1321,9 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
...
@@ -1288,11 +1321,9 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
double
a
=
norm
(
va
),
b
=
norm
(
vb
),
c
=
norm
(
vc
);
double
a
=
norm
(
va
),
b
=
norm
(
vb
),
c
=
norm
(
vc
);
SVector3
n
=
crossprod
(
va
,
vb
);
SVector3
n
=
crossprod
(
va
,
vb
);
//- compute theta
double
theta
=
acos
((
a
*
a
+
b
*
b
-
c
*
c
)
/
(
2
*
a
*
b
));
//in rad
double
theta
=
acos
((
a
*
a
+
b
*
b
-
c
*
c
)
/
(
2
*
a
*
b
));
//in rad
double
sign
=
1.
;
double
sign
=
1.
;
if
(
n
.
z
()
<
0.0
)
{
sign
=
-
1.
;
theta
=
2
*
M_PI
-
theta
;}
if
(
n
.
z
()
<
0.0
)
{
sign
=
-
1.
;
theta
=
2
*
M_PI
-
theta
;}
//printf("theta in deg =%g \n", theta*180./M_PI);
sumTheta
+=
theta
;
sumTheta
+=
theta
;
//- compute grad_uv theta
//- compute grad_uv theta
...
@@ -1300,58 +1331,43 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
...
@@ -1300,58 +1331,43 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
//a=sqrt((u2-u1)^2+(v2-v1)^2))
//a=sqrt((u2-u1)^2+(v2-v1)^2))
//b=sqrt((u3-u2)^2+(v3-v2)^2))
//b=sqrt((u3-u2)^2+(v3-v2)^2))
//c=sqrt((u1-u3)^2+(v1-v3)^2))
//c=sqrt((u1-u3)^2+(v1-v3)^2))
double
u1
=
p1
.
x
();
double
v1
=
p1
.
y
();
double
dTdu1
,
dTdv1
,
dTdu2
,
dTdv2
,
dTdu3
,
dTdv3
;
double
u2
=
p2
.
x
();
double
v2
=
p2
.
y
();
computeThetaDerivatives
(
prev
,
curr
,
next
,
double
u3
=
p3
.
x
();
double
v3
=
p3
.
y
();
dTdu1
,
dTdv1
,
dTdu2
,
dTdv2
,
dTdu3
,
dTdv3
);
double
dTdu1
=
(
v1
*
v3
*
u1
-
v1
*
v3
*
u2
-
v3
*
v2
*
u1
+
v3
*
v2
*
u2
-
v2
*
v1
*
u1
+
2
*
u3
*
v2
*
v1
+
u2
*
SQU
(
v1
)
-
u3
*
SQU
(
v2
)
-
u3
*
SQU
(
v1
)
+
SQU
(
v2
)
*
u1
-
v2
*
v1
*
u2
)
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
sqrt
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
))
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
double
dTdv1
=
(
v1
*
u1
*
u3
-
SQU
(
u1
)
*
v3
+
v1
*
SQU
(
u2
)
-
u2
*
u3
*
v1
-
v1
*
u1
*
u2
-
u2
*
u1
*
v2
+
u3
*
u2
*
v2
-
u1
*
u3
*
v2
+
2
*
v3
*
u2
*
u1
-
v3
*
SQU
(
u2
)
+
v2
*
SQU
(
u1
))
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
sqrt
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
))
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
double
dTdu2
=
(
u1
*
SQU
(
v1
)
*
SQU
(
v2
)
+
u1
*
SQU
(
v1
)
*
SQU
(
v3
)
+
v2
*
CUB
(
v1
)
*
u3
-
u2
*
v2
*
CUB
(
v1
)
-
v3
*
CUB
(
v2
)
*
u3
+
CUB
(
u1
)
*
SQU
(
v3
)
+
CUB
(
u1
)
*
SQU
(
v2
)
+
5
*
u2
*
v2
*
v1
*
SQU
(
u3
)
-
3
*
SQU
(
u2
)
*
v2
*
v1
*
u3
+
u2
*
v2
*
v1
*
SQU
(
v3
)
-
2
*
u2
*
SQU
(
v2
)
*
v1
*
v3
+
u1
*
v2
*
v1
*
SQU
(
u3
)
+
3
*
u1
*
v2
*
v1
*
SQU
(
u2
)
+
u1
*
v2
*
v1
*
SQU
(
v3
)
+
u1
*
SQU
(
v2
)
*
v1
*
v3
+
4
*
u3
*
u2
*
u1
*
SQU
(
v2
)
-
u2
*
v3
*
v2
*
SQU
(
u3
)
+
3
*
SQU
(
u2
)
*
v3
*
v2
*
u3
+
u1
*
v3
*
v2
*
SQU
(
u3
)
-
3
*
u1
*
v3
*
v2
*
SQU
(
u2
)
-
u1
*
v1
*
v3
*
SQU
(
u3
)
+
v2
*
v1
*
u3
*
SQU
(
u1
)
-
u2
*
v2
*
v1
*
SQU
(
u1
)
+
v3
*
SQU
(
v2
)
*
u3
*
v1
+
v3
*
v2
*
u3
*
SQU
(
v1
)
-
3
*
u3
*
v1
*
v3
*
SQU
(
u2
)
-
u3
*
v1
*
v3
*
SQU
(
u1
)
+
u2
*
v1
*
v3
*
SQU
(
u1
)
-
4
*
u1
*
v2
*
v1
*
u3
*
u2
-
u1
*
CUB
(
v2
)
*
v1
-
u2
*
CUB
(
v3
)
*
v2
+
u2
*
SQU
(
v3
)
*
SQU
(
v2
)
+
u1
*
CUB
(
v3
)
*
v2
-
2
*
u1
*
SQU
(
v3
)
*
SQU
(
v2
)
+
u1
*
v3
*
CUB
(
v2
)
-
u1
*
v1
*
CUB
(
v3
)
-
SQU
(
v2
)
*
u3
*
SQU
(
u1
)
+
CUB
(
v2
)
*
u3
*
v1
-
2
*
SQU
(
v2
)
*
u3
*
SQU
(
v1
)
+
3
*
u3
*
SQU
(
u2
)
*
SQU
(
v1
)
+
4
*
u1
*
v1
*
v3
*
u3
*
u2
-
3
*
u1
*
v1
*
v3
*
SQU
(
u2
)
+
SQU
(
u1
)
*
u3
*
v3
*
v2
+
u2
*
SQU
(
v1
)
*
v3
*
v2
+
5
*
u2
*
SQU
(
u1
)
*
v3
*
v2
-
2
*
u1
*
SQU
(
v1
)
*
v3
*
v2
-
4
*
u3
*
u2
*
u1
*
v3
*
v2
-
2
*
u2
*
SQU
(
v2
)
*
SQU
(
u1
)
+
u2
*
SQU
(
v2
)
*
SQU
(
v1
)
-
2
*
SQU
(
u3
)
*
u2
*
SQU
(
v2
)
-
3
*
SQU
(
u3
)
*
u2
*
SQU
(
v1
)
-
u1
*
SQU
(
u3
)
*
SQU
(
v2
)
+
3
*
SQU
(
u2
)
*
u1
*
SQU
(
v3
)
+
SQU
(
v1
)
*
CUB
(
u3
)
-
2
*
u2
*
SQU
(
v1
)
*
SQU
(
v3
)
-
3
*
u2
*
SQU
(
u1
)
*
SQU
(
v3
)
+
u2
*
v1
*
v3
*
SQU
(
u3
)
-
2
*
u3
*
v2
*
v1
*
SQU
(
v3
)
-
2
*
CUB
(
u1
)
*
v3
*
v2
-
CUB
(
u2
)
*
SQU
(
v3
)
-
SQU
(
v1
)
*
CUB
(
u2
)
-
u3
*
CUB
(
v1
)
*
v3
+
2
*
CUB
(
u2
)
*
v1
*
v3
+
u2
*
CUB
(
v1
)
*
v3
+
u2
*
v1
*
CUB
(
v3
)
-
2
*
v2
*
v1
*
CUB
(
u3
)
+
u3
*
SQU
(
v1
)
*
SQU
(
v3
)
+
u3
*
SQU
(
v2
)
*
SQU
(
v3
)
+
SQU
(
v2
)
*
CUB
(
u3
))
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
double
dTdv2
=
-
(
CUB
(
u3
)
*
u2
*
v2
-
u1
*
CUB
(
u3
)
*
v2
-
2
*
u1
*
u3
*
CUB
(
v2
)
+
3
*
u1
*
u3
*
SQU
(
v2
)
*
v1
-
u1
*
u3
*
v2
*
SQU
(
v1
)
+
v3
*
u1
*
u3
*
SQU
(
v1
)
-
4
*
v3
*
u1
*
u3
*
v2
*
v1
+
v1
*
u1
*
u3
*
SQU
(
v3
)
+
u3
*
u2
*
v2
*
SQU
(
v3
)
-
u2
*
u3
*
v1
*
SQU
(
v3
)
+
v3
*
CUB
(
u1
)
*
u3
-
u2
*
CUB
(
u3
)
*
v1
+
u2
*
CUB
(
u1
)
*
v2
+
v1
*
u1
*
CUB
(
u3
)
+
v1
*
SQU
(
u2
)
*
SQU
(
v3
)
+
v1
*
CUB
(
u2
)
*
u1
-
v1
*
CUB
(
u2
)
*
u3
+
2
*
v1
*
SQU
(
u2
)
*
SQU
(
u3
)
-
CUB
(
u2
)
*
u1
*
v3
-
v2
*
SQU
(
u1
)
*
u3
*
u2
+
2
*
v2
*
SQU
(
u2
)
*
u1
*
u3
-
v2
*
u2
*
u1
*
SQU
(
u3
)
-
u1
*
u3
*
v2
*
SQU
(
v3
)
+
3
*
u1
*
u3
*
SQU
(
v2
)
*
v3
+
2
*
SQU
(
u2
)
*
SQU
(
u1
)
*
v3
-
u2
*
CUB
(
u1
)
*
v3
+
4
*
u2
*
u1
*
v3
*
v2
*
v1
-
CUB
(
v1
)
*
SQU
(
u3
)
-
CUB
(
v1
)
*
SQU
(
u2
)
+
2
*
CUB
(
v1
)
*
u3
*
u2
+
2
*
v2
*
SQU
(
v1
)
*
SQU
(
u2
)
+
3
*
v2
*
SQU
(
v1
)
*
SQU
(
u3
)
-
3
*
v1
*
SQU
(
v2
)
*
SQU
(
u3
)
-
v1
*
SQU
(
u1
)
*
SQU
(
u2
)
-
v1
*
SQU
(
u1
)
*
SQU
(
u3
)
+
3
*
SQU
(
u1
)
*
SQU
(
v3
)
*
v2
+
2
*
u2
*
u1
*
CUB
(
v3
)
+
2
*
SQU
(
u2
)
*
SQU
(
v3
)
*
v2
+
v3
*
SQU
(
v1
)
*
SQU
(
u2
)
-
3
*
v3
*
SQU
(
u1
)
*
SQU
(
v2
)
-
5
*
u2
*
u1
*
SQU
(
v3
)
*
v2
+
u2
*
u1
*
v2
*
SQU
(
v1
)
-
SQU
(
u1
)
*
CUB
(
v3
)
-
SQU
(
u2
)
*
CUB
(
v3
)
+
4
*
v3
*
v2
*
v1
*
u3
*
u2
-
v3
*
SQU
(
u1
)
*
SQU
(
u3
)
+
v3
*
CUB
(
u2
)
*
u3
-
v3
*
SQU
(
u2
)
*
SQU
(
u3
)
-
SQU
(
u2
)
*
v2
*
SQU
(
u3
)
+
3
*
u2
*
u1
*
v3
*
SQU
(
v2
)
-
u2
*
u1
*
v3
*
SQU
(
v1
)
-
4
*
v3
*
v2
*
v1
*
SQU
(
u2
)
-
5
*
v2
*
SQU
(
v1
)
*
u3
*
u2
+
3
*
v1
*
SQU
(
v2
)
*
u3
*
u2
-
CUB
(
u1
)
*
u3
*
v2
-
v2
*
SQU
(
u1
)
*
SQU
(
u2
)
+
2
*
v2
*
SQU
(
u1
)
*
SQU
(
u3
)
+
SQU
(
u1
)
*
CUB
(
v2
)
+
CUB
(
v2
)
*
SQU
(
u3
)
+
2
*
v1
*
SQU
(
u1
)
*
u3
*
u2
-
3
*
v1
*
u2
*
u1
*
SQU
(
v2
)
-
v1
*
u2
*
u1
*
SQU
(
v3
)
-
v1
*
SQU
(
u2
)
*
u1
*
u3
-
v1
*
u2
*
u1
*
SQU
(
u3
)
-
v3
*
SQU
(
v1
)
*
u3
*
u2
-
3
*
v3
*
SQU
(
v2
)
*
u3
*
u2
-
v3
*
SQU
(
u1
)
*
u3
*
u2
-
v3
*
SQU
(
u2
)
*
u1
*
u3
+
2
*
v3
*
u2
*
u1
*
SQU
(
u3
))
/
pow
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
),
3.
/
2.
)
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
double
dTdu3
=
-
(
-
SQU
(
v2
)
*
u3
-
2
*
u1
*
v3
*
v2
-
u2
*
SQU
(
v3
)
+
u2
*
v3
*
v2
+
v2
*
v1
*
u3
+
u2
*
v1
*
v3
+
u1
*
SQU
(
v3
)
-
u2
*
v2
*
v1
+
u1
*
SQU
(
v2
)
+
v3
*
v2
*
u3
-
u3
*
v1
*
v3
)
/
sqrt
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
double
dTdv3
=
(
v3
*
u1
*
u3
-
v1
*
SQU
(
u3
)
+
2
*
u2
*
u3
*
v1
-
v1
*
SQU
(
u2
)
+
u2
*
u1
*
v2
-
u2
*
u1
*
v3
-
u3
*
u2
*
v2
-
u3
*
u2
*
v3
-
u1
*
u3
*
v2
+
SQU
(
u2
)
*
v3
+
v2
*
SQU
(
u3
))
/
sqrt
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
pow
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
),
3.
/
2.
)
/
sqrt
(
SQU
(
u2
*
v3
-
u1
*
v3
+
u1
*
v2
-
u2
*
v1
-
v2
*
u3
+
v1
*
u3
)
/
(
SQU
(
u2
)
-
2
*
u2
*
u1
+
SQU
(
u1
)
+
SQU
(
v2
)
-
2
*
v2
*
v1
+
SQU
(
v1
))
/
(
SQU
(
u3
)
-
2
*
u3
*
u2
+
SQU
(
u2
)
+
SQU
(
v3
)
-
2
*
v3
*
v2
+
SQU
(
v2
)));
//- assemble constraint terms
//- assemble constraint terms
myAssembler
.
assemble
(
lag
,
0
,
3
,
prev
,
0
,
1
,
-
dTdu1
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
prev
,
0
,
1
,
-
fac
*
dTdu1
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
prev
,
0
,
2
,
-
dTdv1
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
prev
,
0
,
2
,
-
fac
*
dTdv1
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
curr
,
0
,
1
,
-
dTdu2
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
curr
,
0
,
1
,
-
fac
*
dTdu2
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
curr
,
0
,
2
,
-
dTdv2
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
curr
,
0
,
2
,
-
fac
*
dTdv2
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
next
,
0
,
1
,
-
dTdu3
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
next
,
0
,
1
,
-
fac
*
dTdu3
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
next
,
0
,
2
,
-
dTdv3
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
next
,
0
,
2
,
-
fac
*
dTdv3
);
myAssembler
.
assemble
(
prev
,
0
,
1
,
lag
,
0
,
3
,
-
dTdu1
);
myAssembler
.
assemble
(
prev
,
0
,
1
,
lag
,
0
,
3
,
-
fac
*
dTdu1
);
myAssembler
.
assemble
(
prev
,
0
,
2
,
lag
,
0
,
3
,
-
dTdv1
);
myAssembler
.
assemble
(
prev
,
0
,
2
,
lag
,
0
,
3
,
-
fac
*
dTdv1
);
myAssembler
.
assemble
(
curr
,
0
,
1
,
lag
,
0
,
3
,
-
dTdu2
);
myAssembler
.
assemble
(
curr
,
0
,
1
,
lag
,
0
,
3
,
-
fac
*
dTdu2
);
myAssembler
.
assemble
(
curr
,
0
,
2
,
lag
,
0
,
3
,
-
dTdv2
);
myAssembler
.
assemble
(
curr
,
0
,
2
,
lag
,
0
,
3
,
-
fac
*
dTdv2
);
myAssembler
.
assemble
(
next
,
0
,
1
,
lag
,
0
,
3
,
-
dTdu3
);
myAssembler
.
assemble
(
next
,
0
,
1
,
lag
,
0
,
3
,
-
fac
*
dTdu3
);
myAssembler
.
assemble
(
next
,
0
,
2
,
lag
,
0
,
3
,
-
dTdv3
);
myAssembler
.
assemble
(
next
,
0
,
2
,
lag
,
0
,
3
,
-
fac
*
dTdv3
);
myAssembler
.
assemble
(
prev
,
0
,
1
,
lambda
*
dTdu1
);
myAssembler
.
assemble
(
prev
,
0
,
1
,
lambda
*
fac
*
dTdu1
);
myAssembler
.
assemble
(
prev
,
0
,
2
,
lambda
*
dTdv1
);
myAssembler
.
assemble
(
prev
,
0
,
2
,
lambda
*
fac
*
dTdv1
);
myAssembler
.
assemble
(
curr
,
0
,
1
,
lambda
*
dTdu2
);
myAssembler
.
assemble
(
curr
,
0
,
1
,
lambda
*
fac
*
dTdu2
);
myAssembler
.
assemble
(
curr
,
0
,
2
,
lambda
*
dTdv2
);
myAssembler
.
assemble
(
curr
,
0
,
2
,
lambda
*
fac
*
dTdv2
);
myAssembler
.
assemble
(
next
,
0
,
1
,
lambda
*
dTdu3
);
myAssembler
.
assemble
(
next
,
0
,
1
,
lambda
*
fac
*
dTdu3
);
myAssembler
.
assemble
(
next
,
0
,
2
,
lambda
*
dTdv3
);
myAssembler
.
assemble
(
next
,
0
,
2
,
lambda
*
fac
*
dTdv3
);
}
}
//--- compute constraint
//--- compute constraint
double
G
=
fabs
(
sumTheta
-
(
nb
-
2
)
*
M_PI
);
double
G
=
sumTheta
-
(
nb
-
2
)
*
M_PI
;
printf
(
"Sum of angles G = %g
\n
"
,
G
);
printf
(
"**NL** Sum of angles G = %g
\n
"
,
G
);
double
small
=
0.0
;
myAssembler
.
assemble
(
lag
,
0
,
3
,
lag
,
0
,
3
,
0.0
);
myAssembler
.
assemble
(
lag
,
0
,
3
,
lag
,
0
,
3
,
small
);
//change this
myAssembler
.
assemble
(
lag
,
0
,
3
,
fac
*
G
);
//--- assemble rhs of linear system
myAssembler
.
assemble
(
lag
,
0
,
3
,
G
);
//--- solve linear system
//--- solve linear system
Msg
::
Debug
(
"Assembly done"
);
Msg
::
Debug
(
"Assembly done"
);
_
lsys
->
systemSolve
();
lsys
NL
->
systemSolve
();
Msg
::
Debug
(
"System solved"
);
Msg
::
Debug
(
"System solved"
);
//-- update newton
//-- update newton
...
@@ -1365,26 +1381,42 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
...
@@ -1365,26 +1381,42 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
myAssembler
.
getDofValue
(
v
,
0
,
1
,
du
);
myAssembler
.
getDofValue
(
v
,
0
,
1
,
du
);
myAssembler
.
getDofValue
(
v
,
0
,
2
,
dv
);
myAssembler
.
getDofValue
(
v
,
0
,
2
,
dv
);
meandu
+=
du
;
meandv
+=
dv
;
meandu
+=
du
;
meandv
+=
dv
;
SPoint3
old
=
coordinates
[
v
];
std
::
map
<
MVertex
*
,
SPoint3
>::
iterator
itf
=
coordinates
.
find
(
v
);
coordinates
[
v
]
=
SPoint3
(
old
.
x
()
+
du
,
old
.
y
()
+
dv
,
0.0
);
double
oldu
=
itf
->
second
.
x
(),
oldv
=
itf
->
second
.
y
();
if
(
itf
!=
coordinates
.
end
()){
itf
->
second
[
0
]
=
oldu
+
du
;
itf
->
second
[
1
]
=
oldv
+
dv
;
}
}
}
meandu
/=
allNodes
.
size
();
meandu
/=
allNodes
.
size
();
meandv
/=
allNodes
.
size
();
meandv
/=
allNodes
.
size
();
double
dLambda
;
double
dLambda
;
myAssembler
.
getDofValue
(
lag
,
0
,
3
,
dLambda
);
myAssembler
.
getDofValue
(
lag
,
0
,
3
,
dLambda
);
lambda
=
lambda
+
dLambda
;
lambda
=
lambda
+
dLambda
;
printf
(
"System solved: meandu=%g, meandv=%g dLambda=%g lambda=%g
\n
"
,
meandu
,
meandv
,
dLambda
,
lambda
);
lsysNL
->
zeroMatrix
();
//_lsys->clear();
lsysNL
->
zeroRightHandSide
();
_lsys
->
zeroMatrix
();
_lsys
->
zeroRightHandSide
();
//--priting
//--priting
printStuff
(
iNewton
);
printStuff
(
iNewton
);
//-- exit newton criteria
bool
noOverlap
=
checkOverlap
(
ordered
);
printf
(
"**NL** ---- iNewton %d ---
\n
"
,
iNewton
);
printf
(
"**NL** System solved: du=%g, dv=%g dL=%g
\n
"
,
meandu
,
meandv
,
dLambda
);
if
(
noOverlap
)
{
if
(
checkOrientation
(
0
)){
converged
=
true
;
break
;
}
}
}
//end Newton
}
//end Newton
exit
(
1
);
lsysNL
->
clear
();
//exit(1);
return
converged
;
}
}
...
@@ -1477,7 +1509,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
...
@@ -1477,7 +1509,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
//mettre max NC contraintes par bord
//mettre max NC contraintes par bord
int
NB
=
ordered
.
size
();
int
NB
=
ordered
.
size
();
int
NC
=
std
::
min
(
5
0
,
NB
);
int
NC
=
std
::
min
(
7
0
,
NB
);
int
jump
=
(
int
)
NB
/
NC
;
int
jump
=
(
int
)
NB
/
NC
;
int
nbLoop
=
(
int
)
NB
/
jump
;
int
nbLoop
=
(
int
)
NB
/
jump
;
//printf("nb bound nodes=%d jump =%d \n", NB, jump);
//printf("nb bound nodes=%d jump =%d \n", NB, jump);
...
@@ -1546,7 +1578,6 @@ bool GFaceCompound::parametrize_conformal() const
...
@@ -1546,7 +1578,6 @@ bool GFaceCompound::parametrize_conformal() const
myAssembler
.
fixVertex
(
v1
,
0
,
2
,
0.
);
myAssembler
.
fixVertex
(
v1
,
0
,
2
,
0.
);
myAssembler
.
fixVertex
(
v2
,
0
,
1
,
-
1.
);
myAssembler
.
fixVertex
(
v2
,
0
,
1
,
-
1.
);
myAssembler
.
fixVertex
(
v2
,
0
,
2
,
0.
);
myAssembler
.
fixVertex
(
v2
,
0
,
2
,
0.
);
//printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
// MVertex *v2 ;
// MVertex *v2 ;
// double maxSize = 0.0;
// double maxSize = 0.0;
...
...
This diff is collapsed.
Click to expand it.
Geo/GFaceCompound.h
+
5
−
0
View file @
5e667aa6
...
@@ -100,6 +100,11 @@ class GFaceCompound : public GFace {
...
@@ -100,6 +100,11 @@ class GFaceCompound : public GFace {
SBoundingBox3d
boundEdges
(
const
std
::
list
<
GEdge
*
>
&
elist
)
const
;
SBoundingBox3d
boundEdges
(
const
std
::
list
<
GEdge
*
>
&
elist
)
const
;
SOrientedBoundingBox
obb_boundEdges
(
const
std
::
list
<
GEdge
*
>
&
elist
)
const
;
SOrientedBoundingBox
obb_boundEdges
(
const
std
::
list
<
GEdge
*
>
&
elist
)
const
;
void
fillNeumannBCS
()
const
;
void
fillNeumannBCS
()
const
;
/* double sumAngles(std::vector<MVertex*> ordered) const; */
void
computeThetaDerivatives
(
MVertex
*
prev
,
MVertex
*
curr
,
MVertex
*
next
,
double
&
dTdu1
,
double
&
dTdv1
,
double
&
dTdu2
,
double
&
dTdv2
,
double
&
dTdu3
,
double
&
dTdv3
)
const
;
public
:
public
:
GFaceCompound
(
GModel
*
m
,
int
tag
,
std
::
list
<
GFace
*>
&
compound
,
GFaceCompound
(
GModel
*
m
,
int
tag
,
std
::
list
<
GFace
*>
&
compound
,
std
::
list
<
GEdge
*>
&
U0
,
std
::
list
<
GEdge
*>
&
U1
,
std
::
list
<
GEdge
*>
&
U0
,
std
::
list
<
GEdge
*>
&
U1
,
...
...
This diff is collapsed.
Click to expand it.
Solver/eigenSolver.cpp
+
1
−
1
View file @
5e667aa6
...
@@ -53,7 +53,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
...
@@ -53,7 +53,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// set some default options
// set some default options
_try
(
EPSSetDimensions
(
eps
,
numEigenValues
,
PETSC_DECIDE
,
PETSC_DECIDE
));
_try
(
EPSSetDimensions
(
eps
,
numEigenValues
,
PETSC_DECIDE
,
PETSC_DECIDE
));
_try
(
EPSSetTolerances
(
eps
,
1.e-
6
,
20
));
//1.e-
7
50
_try
(
EPSSetTolerances
(
eps
,
1.e-
7
,
20
));
//1.e-
6
50
//_try(EPSSetType(eps, EPSKRYLOVSCHUR)); //default
//_try(EPSSetType(eps, EPSKRYLOVSCHUR)); //default
_try
(
EPSSetType
(
eps
,
EPSARNOLDI
));
_try
(
EPSSetType
(
eps
,
EPSARNOLDI
));
//_try(EPSSetType(eps,EPSARPACK));
//_try(EPSSetType(eps,EPSARPACK));
...
...
This diff is collapsed.
Click to expand it.
benchmarks/stl/artery.geo
+
1
−
3
View file @
5e667aa6
...
@@ -3,10 +3,8 @@ Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
...
@@ -3,10 +3,8 @@ Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
Mesh
.
CharacteristicLengthFactor
=
0.05
;
Mesh
.
CharacteristicLengthFactor
=
0.05
;
//merge the stl triangulation
Merge
"artery.stl"
;
Merge
"artery.stl"
;
CreateTopology
;
CreateTopology
;
Compound
Surface
(
100
)
=
{
1
}
;
Compound
Surface
(
100
)
=
{
1
};
Physical
Surface
(
101
)
=
{
100
};
Physical
Surface
(
101
)
=
{
100
};
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
