49 #include "Panzer_Workset_Builder.hpp" 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> 60 #include "Shards_CellTopology.hpp" 61 #include "Intrepid2_FunctionSpaceTools.hpp" 62 #include "Intrepid2_CellTools.hpp" 63 #include "Teuchos_Assert.hpp" 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> localSideTopoIDs;
95 std::vector<stk::mesh::Entity> parentElements;
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] <<
")" 111 std::unordered_map<unsigned,std::vector<double> > nodeNormals;
113 TEUCHOS_ASSERT(sides.size() == localSideTopoIDs.size());
114 TEUCHOS_ASSERT(localSideTopoIDs.size() == parentElements.size());
116 RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
117 Intrepid2::DefaultCubatureFactory<double> cubFactory;
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) {
125 std::vector<stk::mesh::Entity> elementEntities;
126 elementEntities.push_back(*parentElement);
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);
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));
142 *pout <<
"element normals: " 143 <<
"gid(" << bulkData->identifier(*parentElement) <<
")" 144 <<
", normal(" <<
normal(0,0,0) <<
"," <<
normal(0,0,1) <<
"," <<
normal(0,0,2) <<
")" 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));
162 for (std::unordered_map<
unsigned,std::vector<double> >::const_iterator node = nodeNormals.begin(); node != nodeNormals.end(); ++node) {
164 TEUCHOS_ASSERT( (node->second.size() % parentTopology->getDimension()) == 0);
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) {
172 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
173 sum += (node->second[surface*parentTopology->getDimension() + dim]) *
174 (node->second[surface*parentTopology->getDimension() + dim]);
176 contribArea[surface] = std::sqrt(
sum);
178 totalArea += contribArea[surface];
182 for (std::size_t i = 0; i < contribArea.size(); ++i)
183 contribArea[i] /= totalArea;
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;
195 *pout <<
"surface normal before normalization: " 196 <<
"gid(" << node->first <<
")" 197 <<
", normal(" <<
normals[node->first][0] <<
"," <<
normals[node->first][1] <<
"," <<
normals[node->first][2] <<
")" 202 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
205 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
209 *pout <<
"surface normal after normalization: " 210 <<
"gid(" << node->first <<
")" 212 <<
normals[node->first][0] <<
"," 213 <<
normals[node->first][1] <<
"," 214 <<
normals[node->first][2] <<
")" 223 const Teuchos::RCP<const panzer_stk::STK_Interface>& mesh,
224 const std::string& sidesetName,
225 const std::string& elementBlockName,
231 std::unordered_map<unsigned,std::vector<double> > nodeEntityIdToNormals;
235 RCP<stk::mesh::MetaData> metaData = mesh->getMetaData();
236 RCP<stk::mesh::BulkData> bulkData = mesh->getBulkData();
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);
247 RCP<const shards::CellTopology> parentTopology = mesh->getCellTopology(elementBlockName);
249 std::vector<std::size_t> localSideTopoIDs;
250 std::vector<stk::mesh::Entity> parentElements;
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) {
259 const size_t numNodes = bulkData->num_nodes(*parentElement);
260 stk::mesh::Entity
const* nodeRelations = bulkData->begin_nodes(*parentElement);
262 normals[mesh->elementLocalId(*parentElement)] = Kokkos::DynRankView<double,PHX::Device>(
"normals",numNodes,parentTopology->getDimension());
264 for (
size_t nodeIndex=0; nodeIndex<numNodes; ++nodeIndex) {
265 stk::mesh::Entity node = nodeRelations[nodeIndex];
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];
274 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
275 (
normals[mesh->elementLocalId(*parentElement)])(nodeIndex,dim) = 0.0;
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.