Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
pugs
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
GitLab community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
code
pugs
Commits
4c68f366
Commit
4c68f366
authored
2 years ago
by
Stéphane Del Pino
Browse files
Options
Downloads
Patches
Plain Diff
[ci-skip] Begin boundary condition treatment
parent
dbcc3085
No related branches found
No related tags found
No related merge requests found
Changes
1
Show whitespace changes
Inline
Side-by-side
Showing
1 changed file
src/mesh/MeshSmootherEscobar.cpp
+189
-213
189 additions, 213 deletions
src/mesh/MeshSmootherEscobar.cpp
with
189 additions
and
213 deletions
src/mesh/MeshSmootherEscobar.cpp
+
189
−
213
View file @
4c68f366
...
...
@@ -40,9 +40,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
using
BoundaryConditionList
=
std
::
vector
<
BoundaryCondition
>
;
BoundaryConditionList
m_boundary_condition_list
;
NodeValue
<
const
bool
>
m_is_fixed_node
;
BoundaryConditionList
_getBCList
(
const
MeshType
&
mesh
,
const
std
::
vector
<
std
::
shared_ptr
<
const
IBoundaryConditionDescriptor
>>&
bc_descriptor_list
)
const
std
::
vector
<
std
::
shared_ptr
<
const
IBoundaryConditionDescriptor
>>&
bc_descriptor_list
)
const
{
BoundaryConditionList
bc_list
;
...
...
@@ -77,8 +79,83 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
return
bc_list
;
}
NodeValue
<
const
bool
>
_getFixedNodes
()
const
{
NodeValue
<
bool
>
is_fixed
{
m_given_mesh
.
connectivity
()};
is_fixed
.
fill
(
false
);
for
(
auto
&&
boundary_condition
:
m_boundary_condition_list
)
{
std
::
visit
(
[
&
](
auto
&&
bc
)
{
using
BCType
=
std
::
decay_t
<
decltype
(
bc
)
>
;
if
constexpr
((
std
::
is_same_v
<
BCType
,
FixedBoundaryCondition
>
)
or
((
Dimension
==
2
)
and
std
::
is_same_v
<
BCType
,
AxisBoundaryCondition
>
))
{
const
Array
<
const
NodeId
>&
node_list
=
bc
.
nodeList
();
parallel_for
(
node_list
.
size
(),
PUGS_LAMBDA
(
size_t
i_node
)
{
is_fixed
[
node_list
[
i_node
]]
=
true
;
});
}
},
boundary_condition
);
}
#warning treat the axis line case in dimension 3 before this
synchronize
(
is_fixed
);
NodeValue
<
int
>
bc_number
{
m_given_mesh
.
connectivity
()};
bc_number
.
fill
(
-
1
);
{
int
bc_id
=
0
;
for
(
auto
&&
boundary_condition
:
m_boundary_condition_list
)
{
std
::
visit
(
[
&
](
auto
&&
bc
)
{
using
BCType
=
std
::
decay_t
<
decltype
(
bc
)
>
;
if
constexpr
(
std
::
is_same_v
<
BCType
,
SymmetryBoundaryCondition
>
)
{
const
Array
<
const
NodeId
>&
node_list
=
bc
.
nodeList
();
parallel_for
(
node_list
.
size
(),
PUGS_LAMBDA
(
size_t
i_node
)
{
const
NodeId
node_id
=
node_list
[
i_node
];
if
(
not
is_fixed
[
node_id
])
{
if
(
bc_number
[
node_id
]
<
0
)
{
bc_number
[
node_id
]
=
bc_id
;
}
else
{
bc_number
[
node_id
]
=
-
2
;
// error tag
}
}
});
}
},
boundary_condition
);
bc_id
++
;
}
if
(
min
(
bc_number
)
<
-
1
)
{
throw
NormalError
(
"Smoothing boundary conditions overlap. One must precise fixed points and smoothing axis"
);
}
}
bool
missing_bc
=
false
;
auto
is_boundary
=
m_given_mesh
.
connectivity
().
isBoundaryNode
();
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
[
&
](
NodeId
node_id
)
{
if
(
is_boundary
[
node_id
])
{
if
(
not
is_fixed
[
node_id
]
and
(
bc_number
[
node_id
]
==
-
1
))
{
missing_bc
=
true
;
}
}
});
if
(
parallel
::
allReduceOr
(
missing_bc
))
{
throw
NormalError
(
"Some boundary conditions are missing (boundary nodes without BC)"
);
}
return
is_fixed
;
}
void
_applyBC
(
NodeValue
<
Rd
>&
shift
)
const
_applyBC
(
const
NodeValue
<
const
Rd
>&
old_xr
,
NodeValue
<
Rd
>&
new_xr
)
const
{
for
(
auto
&&
boundary_condition
:
m_boundary_condition_list
)
{
std
::
visit
(
...
...
@@ -96,22 +173,12 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
node_list
.
size
(),
PUGS_LAMBDA
(
const
size_t
i_node
)
{
const
NodeId
node_id
=
node_list
[
i_node
];
shift
[
node_id
]
=
P
*
shift
[
node_id
];
new_xr
[
node_id
]
=
P
*
new_xr
[
node_id
];
});
}
else
if
constexpr
(
std
::
is_same_v
<
BCType
,
AxisBoundaryCondition
>
)
{
if
constexpr
(
Dimension
>
1
)
{
const
Rd
&
t
=
bc
.
direction
();
const
Rdxd
txt
=
tensorProduct
(
t
,
t
);
const
Array
<
const
NodeId
>&
node_list
=
bc
.
nodeList
();
parallel_for
(
node_list
.
size
(),
PUGS_LAMBDA
(
const
size_t
i_node
)
{
const
NodeId
node_id
=
node_list
[
i_node
];
shift
[
node_id
]
=
txt
*
shift
[
node_id
];
});
throw
NotImplementedError
(
"Escobar: axis boundary conditions"
);
}
else
{
throw
UnexpectedError
(
"AxisBoundaryCondition make no sense in dimension 1"
);
}
...
...
@@ -121,7 +188,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
parallel_for
(
node_list
.
size
(),
PUGS_LAMBDA
(
const
size_t
i_node
)
{
const
NodeId
node_id
=
node_list
[
i_node
];
shift
[
node_id
]
=
zero
;
new_xr
[
node_id
]
=
old_xr
[
node_id
]
;
});
}
else
{
...
...
@@ -132,94 +199,32 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
}
}
NodeValue
<
Rd
>
_getDisplacement
()
const
{
const
ConnectivityType
&
connectivity
=
m_given_mesh
.
connectivity
();
NodeValue
<
const
Rd
>
given_xr
=
m_given_mesh
.
xr
();
auto
node_to_cell_matrix
=
connectivity
.
nodeToCellMatrix
();
auto
cell_to_node_matrix
=
connectivity
.
cellToNodeMatrix
();
auto
node_number_in_their_cells
=
connectivity
.
nodeLocalNumbersInTheirCells
();
NodeValue
<
double
>
max_delta_xr
{
connectivity
};
parallel_for
(
connectivity
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
const
Rd
&
x0
=
given_xr
[
node_id
];
const
auto
&
node_cell_list
=
node_to_cell_matrix
[
node_id
];
double
min_distance_2
=
std
::
numeric_limits
<
double
>::
max
();
for
(
size_t
i_cell
=
0
;
i_cell
<
node_cell_list
.
size
();
++
i_cell
)
{
const
size_t
i_cell_node
=
node_number_in_their_cells
(
node_id
,
i_cell
);
const
CellId
cell_id
=
node_cell_list
[
i_cell
];
const
auto
&
cell_node_list
=
cell_to_node_matrix
[
cell_id
];
for
(
size_t
i_node
=
0
;
i_node
<
cell_node_list
.
size
();
++
i_node
)
{
if
(
i_node
!=
i_cell_node
)
{
const
NodeId
cell_node_id
=
cell_node_list
[
i_node
];
const
Rd
delta
=
x0
-
given_xr
[
cell_node_id
];
min_distance_2
=
std
::
min
(
min_distance_2
,
dot
(
delta
,
delta
));
}
}
}
double
max_delta
=
std
::
sqrt
(
min_distance_2
);
max_delta_xr
[
node_id
]
=
max_delta
;
});
NodeValue
<
Rd
>
shift_r
{
connectivity
};
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
const
auto
&
node_cell_list
=
node_to_cell_matrix
[
node_id
];
Rd
mean_position
(
zero
);
size_t
number_of_neighbours
=
0
;
for
(
size_t
i_cell
=
0
;
i_cell
<
node_cell_list
.
size
();
++
i_cell
)
{
const
size_t
i_cell_node
=
node_number_in_their_cells
(
node_id
,
i_cell
);
const
CellId
cell_id
=
node_cell_list
[
i_cell
];
const
auto
&
cell_node_list
=
cell_to_node_matrix
[
cell_id
];
for
(
size_t
i_node
=
0
;
i_node
<
cell_node_list
.
size
();
++
i_node
)
{
if
(
i_node
!=
i_cell_node
)
{
const
NodeId
cell_node_id
=
cell_node_list
[
i_node
];
mean_position
+=
given_xr
[
cell_node_id
];
number_of_neighbours
++
;
}
}
}
mean_position
=
1.
/
number_of_neighbours
*
mean_position
;
shift_r
[
node_id
]
=
mean_position
-
given_xr
[
node_id
];
});
this
->
_applyBC
(
shift_r
);
synchronize
(
shift_r
);
return
shift_r
;
}
public
:
std
::
shared_ptr
<
const
ItemValueVariant
>
getQuality
()
const
{
if
constexpr
(
Dimension
==
2
)
{
return
std
::
make_shared
<
ItemValueVariant
>
(
m_given_mesh
.
xr
());
}
void
_getNewPositions
(
const
NodeValue
<
const
Rd
>&
old_xr
,
NodeValue
<
Rd
>&
new_xr
)
const
{
const
ConnectivityType
&
connectivity
=
m_given_mesh
.
connectivity
();
NodeValue
<
const
Rd
>
xr
=
m_given_mesh
.
xr
();
auto
is_boundary_node
=
connectivity
.
isBoundaryNode
();
auto
node_is_owned
=
connectivity
.
nodeIsOwned
();
if
constexpr
(
Dimension
==
2
)
{
auto
cell_to_node_matrix
=
connectivity
.
cellToNodeMatrix
();
auto
node_to_cell_matrix
=
connectivity
.
nodeToCellMatrix
();
auto
node_number_in_their_cells
=
connectivity
.
nodeLocalNumbersInTheirCells
();
auto
is_boundary_node
=
connectivity
.
isBoundaryNode
();
NodeValue
<
double
>
quality
{
connectivity
};
constexpr
double
eps
=
1E-15
;
quality
.
fill
(
2
);
auto
f_inner
=
[
=
](
const
NodeId
node_id
,
TinyVector
<
Dimension
>&
x
)
->
double
{
auto
smooth
=
[
=
](
const
NodeId
node_id
,
TinyVector
<
Dimension
>&
x
)
->
double
{
auto
cell_list
=
node_to_cell_matrix
[
node_id
];
auto
node_number_in_cell
=
node_number_in_their_cells
[
node_id
];
...
...
@@ -241,8 +246,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
auto
cell_node_list
=
cell_to_node_matrix
[
cell_list
[
i_cell
]];
const
size_t
cell_nb_nodes
=
cell_node_list
.
size
();
const
TinyVector
a
=
xr
[
cell_node_list
[(
i_cell_node
+
1
)
%
cell_nb_nodes
]];
const
TinyVector
b
=
xr
[
cell_node_list
[(
i_cell_node
+
cell_nb_nodes
-
1
)
%
cell_nb_nodes
]];
const
TinyVector
a
=
old_
xr
[
cell_node_list
[(
i_cell_node
+
1
)
%
cell_nb_nodes
]];
const
TinyVector
b
=
old_
xr
[
cell_node_list
[(
i_cell_node
+
cell_nb_nodes
-
1
)
%
cell_nb_nodes
]];
const
TinyMatrix
<
Dimension
>
A
{
a
[
0
]
-
x
[
0
],
b
[
0
]
-
x
[
0
],
//
a
[
1
]
-
x
[
1
],
b
[
1
]
-
x
[
1
]};
...
...
@@ -262,9 +267,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
auto
frobenius
=
[](
const
TinyMatrix
<
Dimension
>&
M
)
{
return
std
::
sqrt
(
trace
(
transpose
(
M
)
*
M
));
};
// TinyVector<Dimension> f_gradient = zero;
// TinyMatrix<Dimension> f_hessian = zero;
double
final_f
=
0
;
for
(
size_t
i_iter
=
0
;
i_iter
<
100
;
++
i_iter
)
{
...
...
@@ -274,8 +276,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
auto
cell_node_list
=
cell_to_node_matrix
[
cell_list
[
i_cell
]];
const
size_t
cell_nb_nodes
=
cell_node_list
.
size
();
const
TinyVector
a
=
xr
[
cell_node_list
[(
i_cell_node
+
1
)
%
cell_nb_nodes
]];
const
TinyVector
b
=
xr
[
cell_node_list
[(
i_cell_node
+
cell_nb_nodes
-
1
)
%
cell_nb_nodes
]];
const
TinyVector
a
=
old_
xr
[
cell_node_list
[(
i_cell_node
+
1
)
%
cell_nb_nodes
]];
const
TinyVector
b
=
old_
xr
[
cell_node_list
[(
i_cell_node
+
cell_nb_nodes
-
1
)
%
cell_nb_nodes
]];
const
TinyMatrix
<
Dimension
>
A
{
a
[
0
]
-
x
[
0
],
b
[
0
]
-
x
[
0
],
//
a
[
1
]
-
x
[
1
],
b
[
1
]
-
x
[
1
]};
...
...
@@ -311,21 +313,10 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
trace
(
Sigma
*
S_gradient
[
1
])};
const
std
::
array
<
TinyMatrix
<
Dimension
>
,
Dimension
>
//
Sigma_gradient_old
{
sigma_gradient
[
0
]
*
inverse
(
S
)
-
inverse
(
S
)
*
S_gradient
[
0
]
*
Sigma
,
sigma_gradient
[
1
]
*
inverse
(
S
)
-
inverse
(
S
)
*
S_gradient
[
1
]
*
Sigma
};
const
std
::
array
<
TinyMatrix
<
Dimension
>
,
Dimension
>
//
Sigma_gradient_new
{
TinyMatrix
<
Dimension
>
{
0
,
1.
/
std
::
sin
(
alpha
-
1.
/
std
::
tan
(
alpha
)),
//
Sigma_gradient
{
TinyMatrix
<
Dimension
>
{
0
,
1.
/
std
::
sin
(
alpha
-
1.
/
std
::
tan
(
alpha
)),
//
0
,
-
1
},
TinyMatrix
<
Dimension
>
{
-
1.
/
std
::
sin
(
alpha
)
+
1.
/
std
::
tan
(
alpha
),
0
,
//
1
,
0
}};
const
auto
Sigma_gradient
=
Sigma_gradient_new
;
std
::
cout
<<
"Sigma_gradient_old[0] = "
<<
Sigma_gradient_old
[
0
]
<<
'\n'
;
std
::
cout
<<
"Sigma_gradient_new[0] = "
<<
Sigma_gradient_new
[
0
]
<<
'\n'
;
std
::
cout
<<
"Sigma_gradient_old[1] = "
<<
Sigma_gradient_old
[
1
]
<<
'\n'
;
std
::
cout
<<
"Sigma_gradient_new[1] = "
<<
Sigma_gradient_new
[
1
]
<<
'\n'
;
// TinyVector<Dimension> h_gradient = h / (h - sigma_list[i_cell]) * sigma_gradient;
TinyVector
<
Dimension
>
g
{
trace
(
transpose
(
S
)
*
S_gradient
[
0
])
/
S_norm2
//
+
trace
(
transpose
(
Sigma
)
*
Sigma_gradient
[
0
])
/
Sigma_norm2
//
...
...
@@ -360,19 +351,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
f_hessian
+=
f_jr_hessian
;
}
std
::
cout
<<
"f = "
<<
f
<<
'\n'
;
std
::
cout
<<
"grad(f) = "
<<
f_gradient
<<
'\n'
;
std
::
cout
<<
"hess(f) = "
<<
f_hessian
<<
" | hess(f)^T -hess(f) = "
<<
transpose
(
f_hessian
)
-
f_hessian
<<
'\n'
;
std
::
cout
<<
"inv(H) = "
<<
inverse
(
f_hessian
)
<<
'\n'
;
std
::
cout
<<
"inv(H)*grad(f) = "
<<
inverse
(
f_hessian
)
*
f_gradient
<<
'\n'
;
std
::
cout
<<
rang
::
fgB
::
yellow
<<
"x = "
<<
x
<<
" -> "
<<
x
-
inverse
(
f_hessian
)
*
f_gradient
<<
rang
::
fg
::
reset
<<
'\n'
;
std
::
cout
<<
rang
::
fgB
::
green
<<
i_iter
<<
": l2Norm(f_gradient) = "
<<
l2Norm
(
f_gradient
)
<<
rang
::
fg
::
reset
<<
'\n'
;
if
(
l2Norm
(
f_gradient
)
<
1E-6
)
{
break
;
}
...
...
@@ -386,132 +364,128 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
parallel_for
(
connectivity
.
numberOfNodes
(),
PUGS_LAMBDA
(
NodeId
node_id
)
{
// auto cell_list = node_to_cell_matrix[node_id];
// auto node_number_in_cell = node_number_in_their_cells[node_id];
if
(
is_boundary_node
[
node_id
])
{
quality
[
node_id
]
=
1
;
}
else
{
TinyVector
x
=
xr
[
node_id
];
quality
[
node_id
]
=
f_inner
(
node_id
,
x
);
std
::
exit
(
0
);
// TinyMatrix<Dimension> B = identity;
if
(
node_is_owned
[
node_id
]
and
not
is_boundary_node
[
node_id
])
{
quality
[
node_id
]
=
smooth
(
node_id
,
new_xr
[
node_id
]);
}
});
return
std
::
make_shared
<
ItemValueVariant
>
(
quality
);
}
else
{
throw
NotImplementedError
(
"Dimension != 2"
);
}
this
->
_applyBC
(
old_xr
,
new_xr
);
}
std
::
shared_ptr
<
const
IMesh
>
getSmoothedMesh
()
const
{
NodeValue
<
const
Rd
>
given_xr
=
m_given_mesh
.
xr
();
NodeValue
<
Rd
>
new_xr
=
copy
(
given_xr
);
NodeValue
<
Rd
>
xr
=
this
->
_getDisplacement
(
);
this
->
_getNewPositions
(
given_xr
,
new_xr
);
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
xr
[
node_id
]
+=
given_xr
[
node_id
];
});
return
std
::
make_shared
<
MeshType
>
(
m_given_mesh
.
shared_connectivity
(),
xr
);
return
std
::
make_shared
<
MeshType
>
(
m_given_mesh
.
shared_connectivity
(),
new_xr
);
}
std
::
shared_ptr
<
const
IMesh
>
getSmoothedMesh
(
const
FunctionSymbolId
&
function_symbol_id
)
const
{
NodeValue
<
const
Rd
>
given_xr
=
m_given_mesh
.
xr
();
//
NodeValue<const Rd> given_xr = m_given_mesh.xr();
NodeValue
<
const
bool
>
is_displaced
=
InterpolateItemValue
<
bool
(
const
Rd
)
>::
interpolate
(
function_symbol_id
,
given_xr
);
//
NodeValue<const bool> is_displaced =
//
InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
NodeValue
<
Rd
>
xr
=
this
->
_getDisplacement
();
//
NodeValue<Rd> xr = this->_getDisplacement();
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
xr
[
node_id
]
=
is_displaced
[
node_id
]
*
xr
[
node_id
]
+
given_xr
[
node_id
];
});
// parallel_for(
// m_given_mesh.numberOfNodes(),
// PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id];
// });
return
std
::
make_shared
<
MeshType
>
(
m_given_mesh
.
shared_connectivity
(),
xr
);
// return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
return
nullptr
;
}
std
::
shared_ptr
<
const
IMesh
>
getSmoothedMesh
(
const
std
::
vector
<
std
::
shared_ptr
<
const
IZoneDescriptor
>>&
zone_descriptor_list
)
const
{
NodeValue
<
const
Rd
>
given_xr
=
m_given_mesh
.
xr
();
auto
node_to_cell_matrix
=
m_given_mesh
.
connectivity
().
nodeToCellMatrix
();
NodeValue
<
bool
>
is_displaced
{
m_given_mesh
.
connectivity
()};
is_displaced
.
fill
(
false
);
for
(
size_t
i_zone
=
0
;
i_zone
<
zone_descriptor_list
.
size
();
++
i_zone
)
{
MeshCellZone
<
Dimension
>
cell_zone
=
getMeshCellZone
(
m_given_mesh
,
*
zone_descriptor_list
[
i_zone
]);
const
auto
cell_list
=
cell_zone
.
cellList
();
CellValue
<
bool
>
is_zone_cell
{
m_given_mesh
.
connectivity
()};
is_zone_cell
.
fill
(
false
);
parallel_for
(
cell_list
.
size
(),
PUGS_LAMBDA
(
const
size_t
i_cell
)
{
is_zone_cell
[
cell_list
[
i_cell
]]
=
true
;
});
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
auto
node_cell_list
=
node_to_cell_matrix
[
node_id
];
bool
displace
=
true
;
for
(
size_t
i_node_cell
=
0
;
i_node_cell
<
node_cell_list
.
size
();
++
i_node_cell
)
{
const
CellId
cell_id
=
node_cell_list
[
i_node_cell
];
displace
&=
is_zone_cell
[
cell_id
];
}
if
(
displace
)
{
is_displaced
[
node_id
]
=
true
;
}
});
}
synchronize
(
is_displaced
);
NodeValue
<
Rd
>
xr
=
this
->
_getDisplacement
();
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
xr
[
node_id
]
=
is_displaced
[
node_id
]
*
xr
[
node_id
]
+
given_xr
[
node_id
];
});
return
std
::
make_shared
<
MeshType
>
(
m_given_mesh
.
shared_connectivity
(),
xr
);
// NodeValue<const Rd> given_xr = m_given_mesh.xr();
// auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
// NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
// is_displaced.fill(false);
// for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
// MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
// const auto cell_list = cell_zone.cellList();
// CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
// is_zone_cell.fill(false);
// parallel_for(
// cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; });
// parallel_for(
// m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
// auto node_cell_list = node_to_cell_matrix[node_id];
// bool displace = true;
// for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
// const CellId cell_id = node_cell_list[i_node_cell];
// displace &= is_zone_cell[cell_id];
// }
// if (displace) {
// is_displaced[node_id] = true;
// }
// });
// }
// synchronize(is_displaced);
// NodeValue<Rd> xr = this->_getDisplacement();
// parallel_for(
// m_given_mesh.numberOfNodes(),
// PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id];
// });
// return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
return
nullptr
;
}
std
::
shared_ptr
<
const
IMesh
>
getSmoothedMesh
(
const
std
::
vector
<
std
::
shared_ptr
<
const
DiscreteFunctionVariant
>>&
discrete_function_variant_list
)
const
{
NodeValue
<
const
Rd
>
given_xr
=
m_given_mesh
.
xr
();
//
NodeValue<const Rd> given_xr = m_given_mesh.xr();
auto
node_to_cell_matrix
=
m_given_mesh
.
connectivity
().
nodeToCellMatrix
();
//
auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
NodeValue
<
bool
>
is_displaced
{
m_given_mesh
.
connectivity
()};
is_displaced
.
fill
(
false
);
//
NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
//
is_displaced.fill(false);
for
(
size_t
i_zone
=
0
;
i_zone
<
discrete_function_variant_list
.
size
();
++
i_zone
)
{
auto
is_zone_cell
=
discrete_function_variant_list
[
i_zone
]
->
get
<
DiscreteFunctionP0
<
Dimension
,
const
double
>>
();
// for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
// auto is_zone_cell = discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const
// double>>();
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
auto
node_cell_list
=
node_to_cell_matrix
[
node_id
];
bool
displace
=
true
;
for
(
size_t
i_node_cell
=
0
;
i_node_cell
<
node_cell_list
.
size
();
++
i_node_cell
)
{
const
CellId
cell_id
=
node_cell_list
[
i_node_cell
];
displace
&=
(
is_zone_cell
[
cell_id
]
!=
0
);
}
if
(
displace
)
{
is_displaced
[
node_id
]
=
true
;
}
});
}
synchronize
(
is_displaced
);
NodeValue
<
Rd
>
xr
=
this
->
_getDisplacement
();
//
parallel_for(
//
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
//
auto node_cell_list = node_to_cell_matrix[node_id];
//
bool displace = true;
//
for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
//
const CellId cell_id = node_cell_list[i_node_cell];
//
displace &= (is_zone_cell[cell_id] != 0);
//
}
//
if (displace) {
//
is_displaced[node_id] = true;
//
}
//
});
//
}
//
synchronize(is_displaced);
//
NodeValue<Rd> xr = this->_getDisplacement();
parallel_for
(
m_given_mesh
.
numberOfNodes
(),
PUGS_LAMBDA
(
const
NodeId
node_id
)
{
xr
[
node_id
]
=
is_displaced
[
node_id
]
*
xr
[
node_id
]
+
given_xr
[
node_id
];
});
// parallel_for(
// m_given_mesh.numberOfNodes(),
// PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id];
// });
// return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
return
std
::
make_shared
<
MeshType
>
(
m_given_mesh
.
shared_connectivity
(),
xr
)
;
return
nullptr
;
}
MeshSmootherEscobar
(
const
MeshSmootherEscobar
&
)
=
delete
;
...
...
@@ -519,7 +493,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
MeshSmootherEscobar
(
const
MeshType
&
given_mesh
,
const
std
::
vector
<
std
::
shared_ptr
<
const
IBoundaryConditionDescriptor
>>&
bc_descriptor_list
)
:
m_given_mesh
(
given_mesh
),
m_boundary_condition_list
(
this
->
_getBCList
(
given_mesh
,
bc_descriptor_list
))
:
m_given_mesh
{
given_mesh
},
m_boundary_condition_list
{
this
->
_getBCList
(
given_mesh
,
bc_descriptor_list
)},
m_is_fixed_node
{
this
->
_getFixedNodes
()}
{}
~
MeshSmootherEscobar
()
=
default
;
...
...
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