17 #include "Panzer_Workset_Builder.hpp"
23 #include <stk_mesh/base/Selector.hpp>
24 #include <stk_mesh/base/GetEntities.hpp>
25 #include <stk_mesh/base/CreateAdjacentEntities.hpp>
27 #include "Shards_CellTopology.hpp"
29 #include "Intrepid2_CellTools_Serial.hpp"
30 #include "Teuchos_Assert.hpp"
32 namespace panzer_stk {
36 const std::string& sidesetName,
37 const std::string& elementBlockName,
53 stk::mesh::Part * sidePart = mesh->
getSideset(sidesetName);
55 stk::mesh::Selector sideSelector = *sidePart;
56 stk::mesh::Selector blockSelector = *elmtPart;
57 stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
58 std::vector<stk::mesh::Entity> sides;
59 stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
61 std::vector<std::size_t> missingElementIndices;
62 std::vector<std::size_t> localSideTopoIDs;
63 std::vector<stk::mesh::Entity> parentElements;
67 std::vector<std::size_t>::reverse_iterator index;
68 for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
69 sides.erase(sides.begin()+*index);
73 for (std::size_t i=0; i < localSideTopoIDs.size(); ++i) {
74 *pout <<
"parent element: "
75 <<
" gid(" << bulkData->identifier(parentElements[i]) <<
")"
76 <<
", local_face(" << localSideTopoIDs[i] <<
")"
85 std::unordered_map<unsigned,std::vector<double> > nodeNormals;
94 std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
95 std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
96 std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
100 auto side_parametrization = Intrepid2::RefSubcellParametrization<Kokkos::HostSpace>::get(2,parentTopology->getKey());
101 Kokkos::DynRankView<double,Kokkos::HostSpace> normal_at_point(
"normal",3);
102 for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
104 std::vector<stk::mesh::Entity> elementEntities;
105 elementEntities.push_back(*parentElement);
108 auto node_view = nodes.get_view();
120 auto jac_at_point = Kokkos::subview(iv.
jac.get_view(), 0, 0, Kokkos::ALL(), Kokkos::ALL());
121 auto jac_at_point_h = Kokkos::create_mirror_view(jac_at_point);
122 Kokkos::deep_copy(jac_at_point_h, jac_at_point);
124 CellTools::Serial::getPhysicalSideNormal(normal_at_point, side_parametrization, jac_at_point_h, *sideID);
128 *pout <<
"element normals: "
129 <<
"gid(" << bulkData->identifier(*parentElement) <<
")"
130 <<
", normal(" << normal_at_point(0) <<
"," << normal_at_point(1) <<
"," << normal_at_point(2) <<
")"
135 const size_t numNodes = bulkData->num_nodes(*side);
136 stk::mesh::Entity
const* nodeRelations = bulkData->begin_nodes(*side);
137 for (
size_t n=0; n<numNodes; ++n) {
138 stk::mesh::Entity node = nodeRelations[n];
139 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
140 nodeNormals[bulkData->identifier(node)].push_back(normal_at_point(dim));
148 for (std::unordered_map<
unsigned,std::vector<double> >::const_iterator node = nodeNormals.begin(); node != nodeNormals.end(); ++node) {
150 TEUCHOS_ASSERT( (node->second.size() % parentTopology->getDimension()) == 0);
152 int numSurfaceContributions = node->second.size() / parentTopology->getDimension();
153 std::vector<double> contribArea(numSurfaceContributions);
154 double totalArea = 0.0;
155 for (
int surface = 0; surface < numSurfaceContributions; ++surface) {
158 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
159 sum += (node->second[surface*parentTopology->getDimension() + dim]) *
160 (node->second[surface*parentTopology->getDimension() + dim]);
162 contribArea[surface] = std::sqrt(sum);
164 totalArea += contribArea[surface];
168 for (std::size_t i = 0; i < contribArea.size(); ++i)
169 contribArea[i] /= totalArea;
172 normals[node->first].resize(parentTopology->getDimension());
173 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
174 normals[node->first][dim] = 0.0;
175 for (
int surface = 0; surface < numSurfaceContributions; ++surface) {
176 normals[node->first][dim] += node->second[surface*parentTopology->getDimension() + dim] * contribArea[surface] / totalArea;
181 *pout <<
"surface normal before normalization: "
182 <<
"gid(" << node->first <<
")"
183 <<
", normal(" << normals[node->first][0] <<
"," << normals[node->first][1] <<
"," << normals[node->first][2] <<
")"
188 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
189 sum += normals[node->first][dim] * normals[node->first][dim];
191 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim)
192 normals[node->first][dim] /= std::sqrt(sum);
195 *pout <<
"surface normal after normalization: "
196 <<
"gid(" << node->first <<
")"
198 << normals[node->first][0] <<
","
199 << normals[node->first][1] <<
","
200 << normals[node->first][2] <<
")"
210 const std::string& sidesetName,
211 const std::string& elementBlockName,
217 std::unordered_map<unsigned,std::vector<double> > nodeEntityIdToNormals;
225 stk::mesh::Part * sidePart = mesh->
getSideset(sidesetName);
227 stk::mesh::Selector sideSelector = *sidePart;
228 stk::mesh::Selector blockSelector = *elmtPart;
229 stk::mesh::Selector mySelector = metaData->universal_part() & blockSelector & sideSelector;
230 std::vector<stk::mesh::Entity> sides;
231 stk::mesh::get_selected_entities(mySelector,bulkData->buckets(metaData->side_rank()),sides);
235 std::vector<std::size_t> missingElementIndices;
236 std::vector<std::size_t> localSideTopoIDs;
237 std::vector<stk::mesh::Entity> parentElements;
241 std::vector<std::size_t>::reverse_iterator index;
242 for(index=missingElementIndices.rbegin(); index != missingElementIndices.rend(); ++index) {
243 sides.erase(sides.begin()+*index);
246 std::vector<stk::mesh::Entity>::const_iterator side = sides.begin();
247 std::vector<std::size_t>::const_iterator sideID = localSideTopoIDs.begin();
248 std::vector<stk::mesh::Entity>::const_iterator parentElement = parentElements.begin();
249 for ( ; sideID != localSideTopoIDs.end(); ++side,++sideID,++parentElement) {
252 const size_t numNodes = bulkData->num_nodes(*parentElement);
253 stk::mesh::Entity
const* nodeRelations = bulkData->begin_nodes(*parentElement);
255 normals[mesh->
elementLocalId(*parentElement)] = Kokkos::DynRankView<double,PHX::Device>(
"normals",numNodes,parentTopology->getDimension());
256 auto normals_h = Kokkos::create_mirror_view(normals[mesh->
elementLocalId(*parentElement)]);
257 for (
size_t nodeIndex=0; nodeIndex<numNodes; ++nodeIndex) {
258 stk::mesh::Entity node = nodeRelations[nodeIndex];
261 if (nodeEntityIdToNormals.find(bulkData->identifier(node)) != nodeEntityIdToNormals.end()) {
262 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
263 normals_h(nodeIndex,dim) = (nodeEntityIdToNormals[bulkData->identifier(node)])[dim];
267 for (
unsigned dim = 0; dim < parentTopology->getDimension(); ++dim) {
268 normals_h(nodeIndex,dim) = 0.0;
272 Kokkos::deep_copy(normals[mesh->
elementLocalId(*parentElement)], normals_h);
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.
void getElementNodesNoResize(const std::vector< std::size_t > &localIds, ArrayT &nodes) const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
stk::mesh::Part * getSideset(const std::string &name) const
std::size_t elementLocalId(stk::mesh::Entity elmt) const
unsigned getDimension() const
get the dimension
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Data for determining cell topology and dimensionality.
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)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
void setupArrays(const Teuchos::RCP< const panzer::IntegrationRule > &ir)
Sizes/allocates memory for arrays.
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
#define TEUCHOS_ASSERT(assertion_test)
void evaluateValues(const PHX::MDField< Scalar, Cell, NODE, Dim > &cell_node_coordinates, const int num_cells=-1, const Teuchos::RCP< const SubcellConnectivity > &face_connectivity=Teuchos::null, const int num_virtual_cells=-1)
Evaluate basis values.
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const