Panzer  Version of the Day
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 
45 #include "Panzer_STK_Interface.hpp"
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.hpp"
63 #include "Teuchos_Assert.hpp"
64 
65 namespace 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 
80  panzer::MDFieldArrayFactory af("",true);
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> localSideTopoIDs;
95  std::vector<stk::mesh::Entity> parentElements;
96  panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements);
97 
98  if (pout != NULL) {
99  for (std::size_t i=0; i < localSideTopoIDs.size(); ++i) {
100  *pout << "parent element: "
101  << " gid(" << bulkData->identifier(parentElements[i]) << ")"
102  << ", local_face(" << localSideTopoIDs[i] << ")"
103  << std::endl;
104  }
105  }
106 
107  // Do a single element at a time so that we don't have to batch up
108  // into faces
109 
110  // maps a panzer local element id to a list of normals
111  std::unordered_map<unsigned,std::vector<double> > nodeNormals;
112 
113  TEUCHOS_ASSERT(sides.size() == localSideTopoIDs.size());
114  TEUCHOS_ASSERT(localSideTopoIDs.size() == parentElements.size());
115 
116  RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
117  Intrepid2::DefaultCubatureFactory<double> cubFactory;
118  int cubDegree = 1;
119 
120  std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
121  std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
122  std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
123  for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
124 
125  std::vector<stk::mesh::Entity> elementEntities;
126  elementEntities.push_back(*parentElement); // notice this is size 1!
127  PHX::MDField<double,panzer::Cell,panzer::NODE,panzer::Dim> vertices
128  = af.buildStaticArray<double,Cell,NODE,Dim>("",elementEntities.size(), parentTopology->getVertexCount(), mesh->getDimension());
129  mesh->getElementVerticesNoResize(elementEntities,elementBlockName,vertices);
130 
131  panzer::CellData sideCellData(1,*sideID,parentTopology); // this is size 1 because elementEntties is size 1!
132  RCP<panzer::IntegrationRule> ir = Teuchos::rcp(new panzer::IntegrationRule(cubDegree,sideCellData));
133 
135  iv.setupArrays(ir);
136  iv.evaluateValues(vertices);
137 
138  Kokkos::DynRankView<double,PHX::Device> normal("normal",1,ir->num_points,parentTopology->getDimension());
139  Intrepid2::CellTools<double>::getPhysicalSideNormals(normal, iv.jac, *sideID, *(ir->topology));
140 
141  if (pout != NULL) {
142  *pout << "element normals: "
143  << "gid(" << bulkData->identifier(*parentElement) << ")"
144  << ", normal(" << normal(0,0,0) << "," << normal(0,0,1) << "," << normal(0,0,2) << ")"
145  << std::endl;
146  }
147 
148  // loop over nodes in nodes in side and add normal contribution for averaging
149  const size_t numNodes = bulkData->num_nodes(*side);
150  stk::mesh::Entity const* nodeRelations = bulkData->begin_nodes(*side);
151  for (size_t n=0; n<numNodes; ++n) {
152  stk::mesh::Entity node = nodeRelations[n];
153  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
154  nodeNormals[bulkData->identifier(node)].push_back(normal(0,0,dim));
155  }
156  }
157 
158  }
159 
160  // Now do the averaging of contributions
161  //std::unordered_map<unsigned,std::vector<double> > normals;
162  for (std::unordered_map<unsigned,std::vector<double> >::const_iterator node = nodeNormals.begin(); node != nodeNormals.end(); ++node) {
163 
164  TEUCHOS_ASSERT( (node->second.size() % parentTopology->getDimension()) == 0);
165 
166  int numSurfaceContributions = node->second.size() / parentTopology->getDimension();
167  std::vector<double> contribArea(numSurfaceContributions);
168  double totalArea = 0.0;
169  for (int surface = 0; surface < numSurfaceContributions; ++surface) {
170 
171  double sum = 0.0;
172  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
173  sum += (node->second[surface*parentTopology->getDimension() + dim]) *
174  (node->second[surface*parentTopology->getDimension() + dim]);
175 
176  contribArea[surface] = std::sqrt(sum);
177 
178  totalArea += contribArea[surface];
179  }
180 
181  // change the contribArea to the scale factor for each contribution
182  for (std::size_t i = 0; i < contribArea.size(); ++i)
183  contribArea[i] /= totalArea;
184 
185  // loop over contributions and compute final normal values
186  normals[node->first].resize(parentTopology->getDimension());
187  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
188  normals[node->first][dim] = 0.0;
189  for (int surface = 0; surface < numSurfaceContributions; ++surface) {
190  normals[node->first][dim] += node->second[surface*parentTopology->getDimension() + dim] * contribArea[surface] / totalArea;
191  }
192  }
193 
194  if (pout != NULL) {
195  *pout << "surface normal before normalization: "
196  << "gid(" << node->first << ")"
197  << ", normal(" << normals[node->first][0] << "," << normals[node->first][1] << "," << normals[node->first][2] << ")"
198  << std::endl;
199  }
200 
201  double sum = 0.0;
202  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
203  sum += normals[node->first][dim] * normals[node->first][dim];
204 
205  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
206  normals[node->first][dim] /= std::sqrt(sum);
207 
208  if (pout != NULL) {
209  *pout << "surface normal after normalization: "
210  << "gid(" << node->first << ")"
211  << ", normal("
212  << normals[node->first][0] << ","
213  << normals[node->first][1] << ","
214  << normals[node->first][2] << ")"
215  << std::endl;
216  }
217 
218  }
219 
220  }
221 
222  void computeSidesetNodeNormals(std::unordered_map<std::size_t,Kokkos::DynRankView<double,PHX::Device> >& normals,
223  const Teuchos::RCP<const panzer_stk::STK_Interface>& mesh,
224  const std::string& sidesetName,
225  const std::string& elementBlockName,
226  std::ostream* out,
227  std::ostream* pout)
228  {
229  using Teuchos::RCP;
230 
231  std::unordered_map<unsigned,std::vector<double> > nodeEntityIdToNormals;
232 
233  computeSidesetNodeNormals(nodeEntityIdToNormals,mesh,sidesetName,elementBlockName,out,pout);
234 
235  RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
236  RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
237 
238  // Grab all nodes for a surface including ghosted to get correct contributions to normal average
239  stk::mesh::Part * sidePart = mesh->getSideset(sidesetName);
240  stk::mesh::Part * elmtPart = mesh->getElementBlockPart(elementBlockName);
241  stk::mesh::Selector sideSelector = *sidePart;
242  stk::mesh::Selector blockSelector = *elmtPart;
243  stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
244  std::vector<stk::mesh::Entity> sides;
245  stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
246 
247  RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
248 
249  std::vector<std::size_t> localSideTopoIDs;
250  std::vector<stk::mesh::Entity> parentElements;
251  panzer_stk::workset_utils::getUniversalSubcellElements(*mesh,elementBlockName,sides,localSideTopoIDs,parentElements);
252 
253  std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
254  std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
255  std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
256  for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
257 
258  // loop over nodes in nodes in side element
259  const size_t numNodes = bulkData->num_nodes(*parentElement);
260  stk::mesh::Entity const* nodeRelations = bulkData->begin_nodes(*parentElement);
261 
262  normals[mesh->elementLocalId(*parentElement)] = Kokkos::DynRankView<double,PHX::Device>("normals",numNodes,parentTopology->getDimension());
263 
264  for (size_t nodeIndex=0; nodeIndex<numNodes; ++nodeIndex) {
265  stk::mesh::Entity node = nodeRelations[nodeIndex];
266  // if the node is on the sideset, insert, otherwise set normal
267  // to zero (it is an interior node of the parent element).
268  if (nodeEntityIdToNormals.find(bulkData->identifier(node)) != nodeEntityIdToNormals.end()) {
269  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
270  (normals[mesh->elementLocalId(*parentElement)])(nodeIndex,dim) = (nodeEntityIdToNormals[bulkData->identifier(node)])[dim];
271  }
272  }
273  else {
274  for (unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
275  (normals[mesh->elementLocalId(*parentElement)])(nodeIndex,dim) = 0.0;
276  }
277  }
278  }
279 
280  }
281 
282  }
283 
284 }
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 *out, std::ostream *pout)
Computes the normals for all nodes associated with a sideset surface.
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &vertex_coordinates)
Cell vertex coordinates, not basis coordinates.
PHX::MDField< ScalarT > sum
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
PHX::MDField< ScalarT > normal
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)
Data for determining cell topology and dimensionality.
PHX::MDField< ScalarT, Cell, Point, Dim > normals
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.