Skip to content

Commit

Permalink
Merge pull request #3850 from roystgnr/constraints_work
Browse files Browse the repository at this point in the history 
Constraints work
  • Loading branch information
@roystgnr
roystgnr committed May 13, 2024
2 parents ba4fcc1 + 0a52283 commit 08baa8f
Show file tree
Hide file tree
Showing 3 changed files with 42 additions and 40 deletions.
8 changes: 4 additions & 4 deletions include/mesh/mesh_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2000,11 +2000,11 @@ class MeshBase : public ParallelObject
// mesh, such as FE nodes whose Rational Bernstein values need to be
// constrained in terms of values on spline control nodes.
//
// _constraint_rows[constrained_node_id][i].first.first is an
// element,
// _constraint_rows[constrained_node_id][i].first.second is the
// _constraint_rows[constrained_node][i].first.first is an
// element (e.g. a NodeElem for a spline control node),
// _constraint_rows[constrained_node][i].first.second is the
// local node id of that element which is a constraining node,
// _constraint_rows[constrained_node_id][i].second is that node's
// _constraint_rows[constrained_node][i].second is that node's
// constraint coefficient.
constraint_rows_type _constraint_rows;

Expand Down
72 changes: 36 additions & 36 deletions src/mesh/exodusII_io.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -549,28 +549,42 @@ void ExodusII_IO::read (const std::string & fname)

}

// The tailing entries in each element's connectivity
// vector are indices to Exodus constraint coefficient
// rows.

// Concatenating these rows gives a matrix with
// all the constraints for the element nodes: each
// column of that matrix is the constraint coefficients
// for the node associated with that column (via the
// Exodus numbering, not the libMesh numbering).
const auto & my_constraint_rows = exio_helper->bex_cv_conn[elem_num];

for (auto elem_node_index :
make_range(elem->n_nodes()))
{
const auto & my_constraint_rows = exio_helper->bex_cv_conn[elem_num];

// New finite element node data: dot product of
// New finite element node data = dot product of
// constraint matrix columns with spline node data.
// Store that column as a key.
// Store each column's non-zero entries, along with
// the global spline node indices, as a key to
// identify shared finite element nodes.
std::vector<std::pair<dof_id_type, Real>> key;

for (auto spline_node_index :
make_range(exio_helper->bex_num_elem_cvs))
{
// Pick out a row of the element constraint matrix
const unsigned long elem_coef_vec_index =
my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based

auto & coef_vec =
bex_constraint_vec(elem_coef_vec_index, *exio_helper);

// Get coef from this node's column intersect that row
const Real coef =
libmesh_vector_at(coef_vec, elem_node_index);

// Get the libMesh node corresponding to that row
const int gi = (elem_num)*exio_helper->bex_num_elem_cvs +
spline_node_index;
const dof_id_type libmesh_node_id =
Expand All @@ -580,47 +594,31 @@ void ExodusII_IO::read (const std::string & fname)
key.emplace_back(libmesh_node_id, coef);
}

// Have we already created this node? Connect it.
if (const auto local_node_it = local_nodes.find(key);
local_node_it != local_nodes.end())
elem->set_node(dyna_elem_defn.nodes[elem_node_index]) = local_node_it->second;
// Have we not yet created this node? Construct it,
// along with its weight and libMesh constraint row,
// then connect it.
else
{
Point p(0);
Real w = 0;
std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>> constraint_row;

for (auto spline_node_index :
make_range(exio_helper->bex_num_elem_cvs))
for (auto [libmesh_spline_node_id, coef] : key)
{
// global => libMesh index, with crazy 1-based data - see comments above
const int gi = (elem_num)*exio_helper->bex_num_elem_cvs +
spline_node_index;
const dof_id_type libmesh_node_id =
exio_helper->node_num_map[exio_helper->connect[gi] - 1] - 1;

const unsigned long elem_coef_vec_index =
my_constraint_rows[spline_node_index] - 1; // Exodus isn't 0-based

const Node & spline_node = mesh.node_ref(libmesh_node_id);

auto & coef_vec =
bex_constraint_vec(elem_coef_vec_index, *exio_helper);
const Real coef =
libmesh_vector_at(coef_vec, elem_node_index);

// We don't need to store 0 entries;
// constraint_rows is a sparse structure.
if (coef)
{
p.add_scaled(spline_node, coef);
const Real spline_w = weights_exist ?
spline_node.get_extra_datum<Real>(weight_index) : 1;
w += coef * spline_w;

const Elem * nodeelem =
libmesh_map_find(spline_nodeelem_ptrs, &spline_node);
constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
}
const Node & spline_node = mesh.node_ref(libmesh_spline_node_id);

p.add_scaled(spline_node, coef);
const Real spline_w = weights_exist ?
spline_node.get_extra_datum<Real>(weight_index) : 1;
w += coef * spline_w;

const Elem * nodeelem =
libmesh_map_find(spline_nodeelem_ptrs, &spline_node);
constraint_row.emplace_back(std::make_pair(nodeelem, 0), coef);
}

Node *n = mesh.add_point(p, n_nodes++);
Expand All @@ -629,7 +627,9 @@ void ExodusII_IO::read (const std::string & fname)

// If we're building disconnected Bezier
// extraction elements then we don't want to
// find the new nodes to reuse later.
// find the new nodes to reuse later; each
// finite element node will connect to only one
// element.
if (!_disc_bex)
local_nodes[key] = n;
elem->set_node(dyna_elem_defn.nodes[elem_node_index]) = n;
Expand Down
2 changes: 2 additions & 0 deletions src/systems/system.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -485,6 +485,8 @@ void System::reinit_constraints()
get_dof_map().create_dof_constraints(_mesh, this->time);
user_constrain();
get_dof_map().process_constraints(_mesh);
if (libMesh::on_command_line ("--print-constraints"))
get_dof_map().print_dof_constraints(libMesh::out);
#endif
get_dof_map().prepare_send_list();
}
Expand Down

0 comments on commit 08baa8f

Please sign in to comment.