68 const Teuchos::RCP<const panzer_stk::STK_Interface>& mesh,
69 const std::string& sidesetName,
70 const std::string& elementBlockName,
82 RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
83 RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
86 stk::mesh::Part * sidePart = mesh->getSideset(sidesetName);
87 stk::mesh::Part * elmtPart = mesh->getElementBlockPart(elementBlockName);
88 stk::mesh::Selector sideSelector = *sidePart;
89 stk::mesh::Selector blockSelector = *elmtPart;
90 stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
91 std::vector<stk::mesh::Entity> sides;
92 stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
94 std::vector<std::size_t> missingElementIndices;
95 std::vector<std::size_t> localSideTopoIDs;
96 std::vector<stk::mesh::Entity> parentElements;
100 std::vector<std::size_t>::reverse_iterator index;
101 for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
102 sides.erase(sides.begin()+*index);
106 for (std::size_t i=0; i < localSideTopoIDs.size(); ++i) {
107 *pout <<
"parent element: "
108 <<
" gid(" << bulkData->identifier(parentElements[i]) <<
")"
109 <<
", local_face(" << localSideTopoIDs[i] <<
")"
118 std::unordered_map<unsigned,std::vector<double> > nodeNormals;
120 TEUCHOS_ASSERT(sides.size() == localSideTopoIDs.size());
121 TEUCHOS_ASSERT(localSideTopoIDs.size() == parentElements.size());
123 RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
127 std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
128 std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
129 std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
133 auto side_parametrization = Intrepid2::RefSubcellParametrization<Kokkos::HostSpace>::get(2,parentTopology->getKey());
134 Kokkos::DynRankView<double,Kokkos::HostSpace> normal_at_point(
"normal",3);
135 for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
137 std::vector<stk::mesh::Entity> elementEntities;
138 elementEntities.push_back(*parentElement);
139 PHX::MDField<double,panzer::Cell,panzer::NODE,panzer::Dim> vertices
140 = af.
buildStaticArray<double,Cell,NODE,Dim>(
"",elementEntities.size(), parentTopology->getVertexCount(), mesh->getDimension());
141 auto vert_view = vertices.get_view();
142 mesh->getElementVerticesNoResize(elementEntities,elementBlockName,vert_view);
153 auto jac_at_point = Kokkos::subview(iv.
jac.get_view(), 0, 0, Kokkos::ALL(), Kokkos::ALL());
154 auto jac_at_point_h = Kokkos::create_mirror_view(jac_at_point);
155 Kokkos::deep_copy(jac_at_point_h, jac_at_point);
157 CellTools::Serial::getPhysicalSideNormal(normal_at_point, side_parametrization, jac_at_point_h, *sideID);
161 *pout <<
"element normals: "
162 <<
"gid(" << bulkData->identifier(*parentElement) <<
")"
163 <<
", normal(" << normal_at_point(0) <<
"," << normal_at_point(1) <<
"," << normal_at_point(2) <<
")"
168 const size_t numNodes = bulkData->num_nodes(*side);
169 stk::mesh::Entity
const* nodeRelations = bulkData->begin_nodes(*side);
170 for (
size_t n=0; n<numNodes; ++n) {
171 stk::mesh::Entity node = nodeRelations[n];
172 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
173 nodeNormals[bulkData->identifier(node)].push_back(normal_at_point(dim));
181 for (std::unordered_map<
unsigned,std::vector<double> >::const_iterator node = nodeNormals.begin(); node != nodeNormals.end(); ++node) {
183 TEUCHOS_ASSERT( (node->second.size() % parentTopology->getDimension()) == 0);
185 int numSurfaceContributions = node->second.size() / parentTopology->getDimension();
186 std::vector<double> contribArea(numSurfaceContributions);
187 double totalArea = 0.0;
188 for (
int surface = 0; surface < numSurfaceContributions; ++surface) {
191 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
192 sum += (node->second[surface*parentTopology->getDimension() + dim]) *
193 (node->second[surface*parentTopology->getDimension() + dim]);
195 contribArea[surface] = std::sqrt(sum);
197 totalArea += contribArea[surface];
201 for (std::size_t i = 0; i < contribArea.size(); ++i)
202 contribArea[i] /= totalArea;
205 normals[node->first].resize(parentTopology->getDimension());
206 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
207 normals[node->first][dim] = 0.0;
208 for (
int surface = 0; surface < numSurfaceContributions; ++surface) {
209 normals[node->first][dim] += node->second[surface*parentTopology->getDimension() + dim] * contribArea[surface] / totalArea;
214 *pout <<
"surface normal before normalization: "
215 <<
"gid(" << node->first <<
")"
216 <<
", normal(" << normals[node->first][0] <<
"," << normals[node->first][1] <<
"," << normals[node->first][2] <<
")"
221 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
222 sum += normals[node->first][dim] * normals[node->first][dim];
224 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
225 normals[node->first][dim] /= std::sqrt(sum);
228 *pout <<
"surface normal after normalization: "
229 <<
"gid(" << node->first <<
")"
231 << normals[node->first][0] <<
","
232 << normals[node->first][1] <<
","
233 << normals[node->first][2] <<
")"
242 const Teuchos::RCP<const panzer_stk::STK_Interface>& mesh,
243 const std::string& sidesetName,
244 const std::string& elementBlockName,
250 std::unordered_map<unsigned,std::vector<double> > nodeEntityIdToNormals;
254 RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
255 RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
258 stk::mesh::Part * sidePart = mesh->getSideset(sidesetName);
259 stk::mesh::Part * elmtPart = mesh->getElementBlockPart(elementBlockName);
260 stk::mesh::Selector sideSelector = *sidePart;
261 stk::mesh::Selector blockSelector = *elmtPart;
262 stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
263 std::vector<stk::mesh::Entity> sides;
264 stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
266 RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
268 std::vector<std::size_t> missingElementIndices;
269 std::vector<std::size_t> localSideTopoIDs;
270 std::vector<stk::mesh::Entity> parentElements;
274 std::vector<std::size_t>::reverse_iterator index;
275 for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
276 sides.erase(sides.begin()+*index);
279 std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
280 std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
281 std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
282 for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
285 const size_t numNodes = bulkData->num_nodes(*parentElement);
286 stk::mesh::Entity
const* nodeRelations = bulkData->begin_nodes(*parentElement);
288 normals[mesh->elementLocalId(*parentElement)] = Kokkos::DynRankView<double,PHX::Device>(
"normals",numNodes,parentTopology->getDimension());
289 auto normals_h = Kokkos::create_mirror_view(normals[mesh->elementLocalId(*parentElement)]);
290 for (
size_t nodeIndex=0; nodeIndex<numNodes; ++nodeIndex) {
291 stk::mesh::Entity node = nodeRelations[nodeIndex];
294 if (nodeEntityIdToNormals.find(bulkData->identifier(node)) != nodeEntityIdToNormals.end()) {
295 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
296 normals_h(nodeIndex,dim) = (nodeEntityIdToNormals[bulkData->identifier(node)])[dim];
300 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
301 normals_h(nodeIndex,dim) = 0.0;
305 Kokkos::deep_copy(normals[mesh->elementLocalId(*parentElement)], normals_h);