Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_SurfaceNodeNormals.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
44
49#include "Panzer_Workset_Builder.hpp"
52#include "Panzer_CellData.hpp"
54
55#include <stk_mesh/base/Selector.hpp>
56#include <stk_mesh/base/GetEntities.hpp>
57#include <stk_mesh/base/GetBuckets.hpp>
58#include <stk_mesh/base/CreateAdjacentEntities.hpp>
59
60#include "Shards_CellTopology.hpp"
61//#include "Intrepid2_FunctionSpaceTools.hpp"
62#include "Intrepid2_CellTools_Serial.hpp"
63#include "Teuchos_Assert.hpp"
64
65namespace panzer_stk {
66
67 void computeSidesetNodeNormals(std::unordered_map<unsigned,std::vector<double> >& normals,
68 const Teuchos::RCP<const panzer_stk::STK_Interface>& mesh,
69 const std::string& sidesetName,
70 const std::string& elementBlockName,
71 std::ostream* /* out */,
72 std::ostream* pout)
73 {
74 using panzer::Cell;
75 using panzer::NODE;
76 using panzer::Dim;
77
78 using Teuchos::RCP;
79
81
82 RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
83 RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
84
85 // Grab all nodes for a surface including ghosted to get correct contributions to normal average
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);
93
94 std::vector<std::size_t> missingElementIndices;
95 std::vector<std::size_t> localSideTopoIDs;
96 std::vector<stk::mesh::Entity> parentElements;
97 panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements,missingElementIndices);
98
99 // Delete all sides whose neighboring element in elementBlockName is not in the current process
100 std::vector<std::size_t>::reverse_iterator index;
101 for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
102 sides.erase(sides.begin()+*index);
103 }
104
105 if (pout != NULL) {
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] << ")"
110 << std::endl;
111 }
112 }
113
114 // Do a single element at a time so that we don't have to batch up
115 // into faces
116
117 // maps a panzer local element id to a list of normals
118 std::unordered_map<unsigned,std::vector<double> > nodeNormals;
119
120 TEUCHOS_ASSERT(sides.size() == localSideTopoIDs.size());
121 TEUCHOS_ASSERT(localSideTopoIDs.size() == parentElements.size());
122
123 RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
124 //Intrepid2::DefaultCubatureFactory cubFactory;
125 int cubDegree = 1;
126
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();
130
131 // KK: invoke serial interface; cubDegree is 1 and integration point is one
132 // for debugging statement, use max dimension
133 auto side_parametrization = Intrepid2::RefSubcellParametrization<Kokkos::HostSpace>::get(2,parentTopology->getKey());
134 Kokkos::DynRankView<double,Kokkos::HostSpace> normal_at_point("normal",3); // parentTopology->getDimension());
135 for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
136
137 std::vector<stk::mesh::Entity> elementEntities;
138 elementEntities.push_back(*parentElement); // notice this is size 1!
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);
143
144 panzer::CellData sideCellData(1,*sideID,parentTopology); // this is size 1 because elementEntties is size 1!
145 RCP<panzer::IntegrationRule> ir = Teuchos::rcp(new panzer::IntegrationRule(cubDegree,sideCellData));
146
148 iv.setupArrays(ir);
149 iv.evaluateValues(vertices);
150
151 // KK: use serial interface; jac_at_point (D,D) from (C,P,D,D)
152 {
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);
156 Intrepid2::Impl::
157 CellTools::Serial::getPhysicalSideNormal(normal_at_point, side_parametrization, jac_at_point_h, *sideID);
158 }
159
160 if (pout != NULL) {
161 *pout << "element normals: "
162 << "gid(" << bulkData->identifier(*parentElement) << ")"
163 << ", normal(" << normal_at_point(0) << "," << normal_at_point(1) << "," << normal_at_point(2) << ")"
164 << std::endl;
165 }
166
167 // loop over nodes in nodes in side and add normal contribution for averaging
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));
174 }
175 }
176
177 }
178
179 // Now do the averaging of contributions
180 //std::unordered_map<unsigned,std::vector<double> > normals;
181 for (std::unordered_map<unsigned,std::vector<double> >::const_iterator node = nodeNormals.begin(); node != nodeNormals.end(); ++node) {
182
183 TEUCHOS_ASSERT( (node->second.size() % parentTopology->getDimension()) == 0);
184
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) {
189
190 double sum = 0.0;
191 for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
192 sum += (node->second[surface*parentTopology->getDimension() + dim]) *
193 (node->second[surface*parentTopology->getDimension() + dim]);
194
195 contribArea[surface] = std::sqrt(sum);
196
197 totalArea += contribArea[surface];
198 }
199
200 // change the contribArea to the scale factor for each contribution
201 for (std::size_t i = 0; i < contribArea.size(); ++i)
202 contribArea[i] /= totalArea;
203
204 // loop over contributions and compute final normal values
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;
210 }
211 }
212
213 if (pout != NULL) {
214 *pout << "surface normal before normalization: "
215 << "gid(" << node->first << ")"
216 << ", normal(" << normals[node->first][0] << "," << normals[node->first][1] << "," << normals[node->first][2] << ")"
217 << std::endl;
218 }
219
220 double sum = 0.0;
221 for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
222 sum += normals[node->first][dim] * normals[node->first][dim];
223
224 for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
225 normals[node->first][dim] /= std::sqrt(sum);
226
227 if (pout != NULL) {
228 *pout << "surface normal after normalization: "
229 << "gid(" << node->first << ")"
230 << ", normal("
231 << normals[node->first][0] << ","
232 << normals[node->first][1] << ","
233 << normals[node->first][2] << ")"
234 << std::endl;
235 }
236
237 }
238
239 }
240
241 void computeSidesetNodeNormals(std::unordered_map<std::size_t,Kokkos::DynRankView<double,PHX::Device> >& normals,
242 const Teuchos::RCP<const panzer_stk::STK_Interface>& mesh,
243 const std::string& sidesetName,
244 const std::string& elementBlockName,
245 std::ostream* out,
246 std::ostream* pout)
247 {
248 using Teuchos::RCP;
249
250 std::unordered_map<unsigned,std::vector<double> > nodeEntityIdToNormals;
251
252 computeSidesetNodeNormals(nodeEntityIdToNormals,mesh,sidesetName,elementBlockName,out,pout);
253
254 RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
255 RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
256
257 // Grab all nodes for a surface including ghosted to get correct contributions to normal average
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);
265
266 RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
267
268 std::vector<std::size_t> missingElementIndices;
269 std::vector<std::size_t> localSideTopoIDs;
270 std::vector<stk::mesh::Entity> parentElements;
271 panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements,missingElementIndices);
272
273 // Delete all sides whose neighboring element in elementBlockName is not in the current process
274 std::vector<std::size_t>::reverse_iterator index;
275 for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
276 sides.erase(sides.begin()+*index);
277 }
278
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) {
283
284 // loop over nodes in nodes in side element
285 const size_t numNodes = bulkData->num_nodes(*parentElement);
286 stk::mesh::Entity const* nodeRelations = bulkData->begin_nodes(*parentElement);
287
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];
292 // if the node is on the sideset, insert, otherwise set normal
293 // to zero (it is an interior node of the parent element).
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];
297 }
298 }
299 else {
300 for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
301 normals_h(nodeIndex,dim) = 0.0;
302 }
303 }
304 }
305 Kokkos::deep_copy(normals[mesh->elementLocalId(*parentElement)], normals_h);
306 }
307
308 }
309
310}
Data for determining cell topology and dimensionality.
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates, const int num_cells=-1, const Teuchos::RCP< const SubcellConnectivity > &face_connectivity=Teuchos::null, const int num_virtual_cells=-1)
Evaluate basis values.
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
void getUniversalSubcellElements(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &entities, std::vector< std::size_t > &localEntityIds, std::vector< stk::mesh::Entity > &elements, std::vector< std::size_t > &missingElementIndices)
void computeSidesetNodeNormals(std::unordered_map< unsigned, std::vector< double > > &normals, const Teuchos::RCP< const panzer_stk::STK_Interface > &mesh, const std::string &sidesetName, const std::string &elementBlockName, std::ostream *, std::ostream *pout)
Computes the normals for all nodes associated with a sideset surface.