On Mon, 27 Feb 2012, Boyce Griffith wrote:
On 2/24/12 7:53 PM, Roy Stogner wrote:
Anything I'm missing? I'm going to start coding and testing
everything this weekend, bearing in mind that I may tear it all up and
start over if Ben's got better ideas based on his recent FIN-S
true-Dirichlet experience.
Let me know if I can help. I have some time this week that I could spend
working on this.
That would be great, since I don't have time to work on this. But
why let that stop me?
Updated patch is posted; thanks also to Derek for offering to test.
Assuming this doesn't break anyone's adaptive, periodic, etc. codes,
then I'll commit and add a Dirichlet constraints generation routine
(probably an API along the lines of the functor version of
System::project_boundary_solution unless anyone dislikes that) this
week.
---
Roy
Index: include/base/dof_map.h
===================================================================
--- include/base/dof_map.h (revision 5180)
+++ include/base/dof_map.h (working copy)
@@ -79,9 +79,9 @@
* a pointer-to-std::map; the destructor isn't virtual!
*/
class DofConstraints : public std::map<unsigned int,
- DofConstraintRow,
+ std::pair<DofConstraintRow,Number>,
std::less<unsigned int>,
- Threads::scalable_allocator<std::pair<const unsigned int, DofConstraintRow> > >
+ Threads::scalable_allocator<std::pair<const unsigned int, std::pair<DofConstraintRow,Number> > > >
{
};
@@ -103,9 +103,9 @@
* a pointer-to-std::map; the destructor isn't virtual!
*/
class NodeConstraints : public std::map<const Node *,
- NodeConstraintRow,
- std::less<const Node *>,
- Threads::scalable_allocator<std::pair<const Node * const, NodeConstraintRow> > >
+ std::pair<NodeConstraintRow,Point>,
+ std::less<const Node *>,
+ Threads::scalable_allocator<std::pair<const Node * const, std::pair<NodeConstraintRow,Point> > > >
{
};
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
@@ -475,12 +475,23 @@
void process_constraints ();
/**
- * Adds a copy of the user-defined row to the constraint matrix.
+ * Adds a copy of the user-defined row to the constraint matrix, using
+ * an inhomogeneous right-hand-side for the constraint equation.
+ */
+ void add_constraint_row (const unsigned int dof_number,
+ const DofConstraintRow& constraint_row,
+ const Number constraint_rhs,
+ const bool forbid_constraint_overwrite);
+
+ /**
+ * Adds a copy of the user-defined row to the constraint matrix, using
+ * a homogeneous right-hand-side for the constraint equation.
* By default, produces an error if the DOF was already constrained.
*/
void add_constraint_row (const unsigned int dof_number,
- const DofConstraintRow& constraint_row,
- const bool forbid_constraint_overwrite = true);
+ const DofConstraintRow& constraint_row,
+ const bool forbid_constraint_overwrite = true)
+ { add_constraint_row(dof_number, constraint_row, 0., forbid_constraint_overwrite); }
/**
* Returns an iterator pointing to the first DoF constraint row
@@ -593,6 +604,28 @@
DenseVector<Number>& rhs,
std::vector<unsigned int>& elem_dofs,
bool asymmetric_constraint_rows = true) const;
+ /**
+ * Constrains the element matrix and vector. This method requires
+ * the element matrix to be square, in which case the elem_dofs
+ * correspond to the global DOF indices of both the rows and
+ * columns of the element matrix. For this case the rows
+ * and columns of the matrix necessarily correspond to variables
+ * of the same approximation order.
+ *
+ * The heterogenous version of this method creates linear systems in
+ * which heterogenously constrained degrees of freedom will solve to
+ * their correct offset values, as would be appropriate for finding
+ * a solution to a linear problem in a single algebraic solve. The
+ * non-heterogenous version of this method creates linear systems in
+ * which even heterogenously constrained degrees of freedom are
+ * solved without offset values taken into account, as would be
+ * appropriate for finding linearized updates to a solution in which
+ * heterogenous constraints are already satisfied.
+ */
+ void heterogenously_constrain_element_matrix_and_vector (DenseMatrix<Number>& matrix,
+ DenseVector<Number>& rhs,
+ std::vector<unsigned int>& elem_dofs,
+ bool asymmetric_constraint_rows = true) const;
/**
* Constrains a dyadic element matrix B = v w'. This method
@@ -823,6 +856,21 @@
const bool called_recursively=false) const;
/**
+ * Build the constraint matrix C and the forcing vector H
+ * associated with the element degree of freedom indices elem_dofs.
+ * The optional parameter \p called_recursively should be left at
+ * the default value \p false. This is used to handle the special
+ * case of an element's degrees of freedom being constrained in
+ * terms of other, local degrees of freedom. The usual case is for
+ * an elements DOFs to be constrained by some other, external DOFs
+ * and/or Dirichlet conditions.
+ */
+ void build_constraint_matrix_and_vector (DenseMatrix<Number>& C,
+ DenseVector<Number>& H,
+ std::vector<unsigned int>& elem_dofs,
+ const bool called_recursively=false) const;
+
+ /**
* Finds all the DOFS associated with the element DOFs elem_dofs.
* This will account for off-element couplings via hanging nodes.
*/
Index: src/fe/fe_base.C
===================================================================
--- src/fe/fe_base.C (revision 5180)
+++ src/fe/fe_base.C (working copy)
@@ -2190,7 +2190,7 @@
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
DofConstraintRow& constraint_row =
- constraints[my_dof_g];
+ constraints[my_dof_g].first;
constraint_row.insert(std::make_pair(their_dof_g,
their_dof_value));
@@ -2325,7 +2325,7 @@
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
// A reference to the constraint row.
- NodeConstraintRow& constraint_row = constraints[my_node];
+ NodeConstraintRow& constraint_row = constraints[my_node].first;
constraint_row.insert(std::make_pair (their_node,
0.));
@@ -2341,7 +2341,7 @@
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
// A reference to the constraint row.
- NodeConstraintRow& constraint_row = constraints[my_node];
+ NodeConstraintRow& constraint_row = constraints[my_node].first;
constraint_row.insert(std::make_pair (their_node,
their_value));
@@ -2571,7 +2571,7 @@
continue;
DofConstraintRow& their_constraint_row =
- constraints[their_dof_g];
+ constraints[their_dof_g].first;
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
@@ -2619,7 +2619,7 @@
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
DofConstraintRow& constraint_row =
- constraints[my_dof_g];
+ constraints[my_dof_g].first;
constraint_row.insert(std::make_pair(their_dof_g,
their_dof_value));
@@ -2777,7 +2777,7 @@
continue;
const NodeConstraintRow& their_constraint_row =
- constraints[their_node];
+ constraints[their_node].first;
for (unsigned int orig_side_n=0;
orig_side_n < n_side_nodes;
@@ -2834,7 +2834,7 @@
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
NodeConstraintRow& constraint_row =
- constraints[my_node];
+ constraints[my_node].first;
constraint_row.insert(std::make_pair(their_node,
their_value));
Index: src/fe/fe_lagrange.C
===================================================================
--- src/fe/fe_lagrange.C (revision 5180)
+++ src/fe/fe_lagrange.C (working copy)
@@ -730,7 +730,7 @@
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
// A reference to the constraint row.
- DofConstraintRow& constraint_row = constraints[my_dof_g];
+ DofConstraintRow& constraint_row = constraints[my_dof_g].first;
constraint_row.insert(std::make_pair (their_dof_g,
their_dof_value));
Index: src/base/dof_map.C
===================================================================
--- src/base/dof_map.C (revision 5180)
+++ src/base/dof_map.C (working copy)
@@ -1973,7 +1973,7 @@
libmesh_assert (pos != _dof_constraints.end());
- const DofConstraintRow& constraint_row = pos->second;
+ const DofConstraintRow& constraint_row = pos->second.first;
// adaptive p refinement currently gives us lots of empty constraint
// rows - we should optimize those DoFs away in the future. [RHS]
@@ -2502,34 +2502,41 @@
#if defined(LIBMESH_ENABLE_AMR) || defined(LIBMESH_ENABLE_PERIODIC)
- unsigned int n_constraints = 0, max_constraint_length = 0;
+ unsigned int n_constraints = 0, max_constraint_length = 0,
+ n_rhss = 0;
long double avg_constraint_length = 0.;
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
- unsigned int n_hanging_nodes = 0, max_node_constraint_length = 0;
+ unsigned int n_node_constraints = 0, max_node_constraint_length = 0,
+ n_node_rhss = 0;
long double avg_node_constraint_length = 0.;
#endif
for (DofConstraints::const_iterator it=_dof_constraints.begin();
it != _dof_constraints.end(); ++it)
{
- const DofConstraintRow& row = it->second;
+ const DofConstraintRow& row = it->second.first;
unsigned int rowsize = row.size();
max_constraint_length = std::max(max_constraint_length,
rowsize);
avg_constraint_length += rowsize;
n_constraints++;
+
+ if (it->second.second)
+ n_rhss++;
}
os << " DofMap Constraints\n Number of DoF Constraints = "
<< n_constraints;
+ if (n_rhss)
+ os << '\n'
+ << " Number of heterogenous constraints= " << n_rhss;
if (n_constraints)
{
avg_constraint_length /= n_constraints;
- os << " DofMap Constraints\n Number of Constraints = " << n_constraints
- << '\n'
+ os << '\n'
<< " Average DoF Constraint Length= " << avg_constraint_length;
}
@@ -2537,19 +2544,25 @@
for (NodeConstraints::const_iterator it=_node_constraints.begin();
it != _node_constraints.end(); ++it)
{
- const NodeConstraintRow& row = it->second;
+ const NodeConstraintRow& row = it->second.first;
unsigned int rowsize = row.size();
max_node_constraint_length = std::max(max_node_constraint_length,
rowsize);
avg_node_constraint_length += rowsize;
- n_hanging_nodes++;
+ n_node_constraints++;
+
+ if (it->second.second != Point(0))
+ n_node_rhss++;
}
- if (n_hanging_nodes)
+ os << "\n Number of Node Constraints = " << n_node_constraints;
+ if (n_rhss)
+ os << '\n'
+ << " Number of heterogenous node constraints= " << n_node_rhss;
+ if (n_node_constraints)
{
- os << "\n Number of Node Constraints = " << n_hanging_nodes;
- avg_node_constraint_length /= n_hanging_nodes;
+ avg_node_constraint_length /= n_node_constraints;
os << "\n Maximum Node Constraint Length= " << max_node_constraint_length
<< '\n'
<< " Average Node Constraint Length= " << avg_node_constraint_length;
Index: src/base/dof_map_constraints.C
===================================================================
--- src/base/dof_map_constraints.C (revision 5180)
+++ src/base/dof_map_constraints.C (working copy)
@@ -276,6 +276,7 @@
void DofMap::add_constraint_row (const unsigned int dof_number,
const DofConstraintRow& constraint_row,
+ const Number constraint_rhs,
const bool forbid_constraint_overwrite)
{
// Optionally allow the user to overwrite constraints. Defaults to false.
@@ -287,7 +288,7 @@
libmesh_error();
}
- std::pair<unsigned int, DofConstraintRow> kv(dof_number, constraint_row);
+ std::pair<unsigned int, std::pair<DofConstraintRow,Number> > kv(dof_number, std::make_pair(constraint_row, constraint_rhs));
_dof_constraints.insert(kv);
}
@@ -304,7 +305,8 @@
it != _node_constraints.end(); ++it)
{
const Node *node = it->first;
- const NodeConstraintRow& row = it->second;
+ const NodeConstraintRow& row = it->second.first;
+ const Point& offset = it->second.second;
os << "Constraints for Node id " << node->id()
<< ": \t";
@@ -314,6 +316,8 @@
os << " (" << pos->first->id() << ","
<< pos->second << ")\t";
+ os << "rhs: " << offset;
+
os << std::endl;
}
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
@@ -325,7 +329,8 @@
it != _dof_constraints.end(); ++it)
{
const unsigned int i = it->first;
- const DofConstraintRow& row = it->second;
+ const DofConstraintRow& row = it->second.first;
+ const Number rhs = it->second.second;
os << "Constraints for DoF " << i
<< ": \t";
@@ -335,6 +340,8 @@
os << " (" << pos->first << ","
<< pos->second << ")\t";
+ os << "rhs: " << rhs;
+
os << std::endl;
}
}
@@ -390,7 +397,7 @@
libmesh_assert (pos != _dof_constraints.end());
- const DofConstraintRow& constraint_row = pos->second;
+ const DofConstraintRow& constraint_row = pos->second.first;
libmesh_assert (!constraint_row.empty());
@@ -463,7 +470,7 @@
libmesh_assert (pos != _dof_constraints.end());
- const DofConstraintRow& constraint_row = pos->second;
+ const DofConstraintRow& constraint_row = pos->second.first;
// p refinement creates empty constraint rows
// libmesh_assert (!constraint_row.empty());
@@ -483,18 +490,97 @@
// compute matrix/vector product
C.vector_mult_transpose(rhs, old_rhs);
+ } // end if is constrained...
- libmesh_assert (elem_dofs.size() == rhs.size());
+ STOP_LOG("cnstrn_elem_mat_vec()", "DofMap");
+}
+
+
+void DofMap::heterogenously_constrain_element_matrix_and_vector
+ (DenseMatrix<Number>& matrix,
+ DenseVector<Number>& rhs,
+ std::vector<unsigned int>& elem_dofs,
+ bool asymmetric_constraint_rows) const
+{
+ libmesh_assert (elem_dofs.size() == matrix.m());
+ libmesh_assert (elem_dofs.size() == matrix.n());
+ libmesh_assert (elem_dofs.size() == rhs.size());
+
+ // check for easy return
+ if (this->n_constrained_dofs() == 0)
+ return;
+
+ // The constrained matrix is built up as C^T K C.
+ // The constrained RHS is built up as C^T (F - K H)
+ DenseMatrix<Number> C;
+ DenseVector<Number> H;
+
+ this->build_constraint_matrix_and_vector (C, H, elem_dofs);
+
+ START_LOG("hetero_cnstrn_elem_mat_vec()", "DofMap");
+
+ // It is possible that the matrix is not constrained at all.
+ if ((C.m() == matrix.m()) &&
+ (C.n() == elem_dofs.size())) // It the matrix is constrained
+ {
+ // Compute matrix/vector product K H
+ DenseVector<Number> KH;
+ matrix.vector_mult(KH, H);
+
+ // Compute the matrix-vector product C^T (F - KH)
+ DenseVector<Number> F_minus_KH(rhs);
+ F_minus_KH -= KH;
+ C.vector_mult_transpose(rhs, F_minus_KH);
+
+ // Compute the matrix-matrix-matrix product C^T K C
+ matrix.left_multiply_transpose (C);
+ matrix.right_multiply (C);
+
+ libmesh_assert (matrix.m() == matrix.n());
+ libmesh_assert (matrix.m() == elem_dofs.size());
+ libmesh_assert (matrix.n() == elem_dofs.size());
+
for (unsigned int i=0; i<elem_dofs.size(); i++)
if (this->is_constrained_dof(elem_dofs[i]))
{
+ for (unsigned int j=0; j<matrix.n(); j++)
+ matrix(i,j) = 0.;
+
// If the DOF is constrained
- rhs(i) = 0.;
+ matrix(i,i) = 1.;
+
+ // This will put a nonsymmetric entry in the constraint
+ // row to ensure that the linear system produces the
+ // correct value for the constrained DOF.
+ if (asymmetric_constraint_rows)
+ {
+ DofConstraints::const_iterator
+ pos = _dof_constraints.find(elem_dofs[i]);
+
+ libmesh_assert (pos != _dof_constraints.end());
+
+ const DofConstraintRow& constraint_row = pos->second.first;
+
+// p refinement creates empty constraint rows
+// libmesh_assert (!constraint_row.empty());
+
+ for (DofConstraintRow::const_iterator
+ it=constraint_row.begin(); it != constraint_row.end();
+ ++it)
+ for (unsigned int j=0; j<elem_dofs.size(); j++)
+ if (elem_dofs[j] == it->first)
+ matrix(i,j) = -it->second;
+
+ rhs(i) = -pos->second.second;
+ }
+ else
+ rhs(i) = 0.;
}
+
} // end if is constrained...
- STOP_LOG("cnstrn_elem_mat_vec()", "DofMap");
+ STOP_LOG("hetero_cnstrn_elem_mat_vec()", "DofMap");
}
@@ -561,7 +647,7 @@
libmesh_assert (pos != _dof_constraints.end());
- const DofConstraintRow& constraint_row = pos->second;
+ const DofConstraintRow& constraint_row = pos->second.first;
libmesh_assert (!constraint_row.empty());
@@ -611,7 +697,12 @@
if (this->is_constrained_dof(row_dofs[i]))
{
// If the DOF is constrained
- rhs(i) = 0.;
+ DofConstraints::const_iterator
+ pos = _dof_constraints.find(row_dofs[i]);
+
+ libmesh_assert (pos != _dof_constraints.end());
+
+ rhs(i) = -pos->second.second;
}
} // end if the RHS is constrained.
@@ -660,7 +751,12 @@
if (this->is_constrained_dof(row_dofs[i]))
{
// If the DOF is constrained
- v(i) = 0.;
+ DofConstraints::const_iterator
+ pos = _dof_constraints.find(row_dofs[i]);
+
+ libmesh_assert (pos != _dof_constraints.end());
+
+ v(i) = -pos->second.second;
}
} // end if the RHS is constrained.
@@ -755,9 +851,9 @@
constrained_dof >= this->end_dof())
continue;
- const DofConstraintRow constraint_row = c_it->second;
+ const DofConstraintRow constraint_row = c_it->second.first;
- Number exact_value = 0;
+ Number exact_value = c_it->second.second;
for (DofConstraintRow::const_iterator
j=constraint_row.begin(); j != constraint_row.end();
++j)
@@ -831,7 +927,12 @@
global_dof >= vec.first_local_index() &&
global_dof < vec.last_local_index())
{
- Number exact_value = 0;
+ DofConstraints::const_iterator
+ pos = _dof_constraints.find(global_dof);
+
+ libmesh_assert (pos != _dof_constraints.end());
+
+ Number exact_value = pos->second.second;
for (unsigned int j=0; j!=C.n(); ++j)
{
if (local_dof_indices[j] != global_dof)
@@ -881,7 +982,7 @@
libmesh_assert (pos != _dof_constraints.end());
- const DofConstraintRow& constraint_row = pos->second;
+ const DofConstraintRow& constraint_row = pos->second.first;
// Constraint rows in p refinement may be empty
// libmesh_assert (!constraint_row.empty());
@@ -933,7 +1034,7 @@
libmesh_assert (pos != _dof_constraints.end());
- const DofConstraintRow& constraint_row = pos->second;
+ const DofConstraintRow& constraint_row = pos->second.first;
// p refinement creates empty constraint rows
// libmesh_assert (!constraint_row.empty());
@@ -972,6 +1073,135 @@
}
+
+void DofMap::build_constraint_matrix_and_vector
+ (DenseMatrix<Number>& C,
+ DenseVector<Number>& H,
+ std::vector<unsigned int>& elem_dofs,
+ const bool called_recursively) const
+{
+ if (!called_recursively)
+ START_LOG("build_constraint_matrix_and_vector()", "DofMap");
+
+ // Create a set containing the DOFs we already depend on
+ typedef std::set<unsigned int> RCSet;
+ RCSet dof_set;
+
+ bool we_have_constraints = false;
+
+ // Next insert any other dofs the current dofs might be constrained
+ // in terms of. Note that in this case we may not be done: Those
+ // may in turn depend on others. So, we need to repeat this process
+ // in that case until the system depends only on unconstrained
+ // degrees of freedom.
+ for (unsigned int i=0; i<elem_dofs.size(); i++)
+ if (this->is_constrained_dof(elem_dofs[i]))
+ {
+ we_have_constraints = true;
+
+ // If the DOF is constrained
+ DofConstraints::const_iterator
+ pos = _dof_constraints.find(elem_dofs[i]);
+
+ libmesh_assert (pos != _dof_constraints.end());
+
+ const DofConstraintRow& constraint_row = pos->second.first;
+
+// Constraint rows in p refinement may be empty
+// libmesh_assert (!constraint_row.empty());
+
+ for (DofConstraintRow::const_iterator
+ it=constraint_row.begin(); it != constraint_row.end();
+ ++it)
+ dof_set.insert (it->first);
+ }
+
+ // May be safe to return at this point
+ // (but remember to stop the perflog)
+ if (!we_have_constraints)
+ {
+ STOP_LOG("build_constraint_matrix_and_vector()", "DofMap");
+ return;
+ }
+
+ // delay inserting elem_dofs for efficiency in the case of
+ // no constraints. In that case we don't get here!
+ dof_set.insert (elem_dofs.begin(),
+ elem_dofs.end());
+
+ // If we added any DOFS then we need to do this recursively.
+ // It is possible that we just added a DOF that is also
+ // constrained!
+ //
+ // Also, we need to handle the special case of an element having DOFs
+ // constrained in terms of other, local DOFs
+ if ((dof_set.size() != elem_dofs.size()) || // case 1: constrained in terms of other DOFs
+ !called_recursively) // case 2: constrained in terms of our own DOFs
+ {
+ // Create a new list of element DOFs containing the
+ // contents of the current dof_set.
+ std::vector<unsigned int> new_elem_dofs (dof_set.begin(),
+ dof_set.end());
+
+ // Now we can build the constraint matrix.
+ // Note that resize also zeros for a DenseMatrix<Number>.
+ C.resize (elem_dofs.size(), new_elem_dofs.size());
+
+ // Create the C constraint matrix.
+ for (unsigned int i=0; i<elem_dofs.size(); i++)
+ if (this->is_constrained_dof(elem_dofs[i]))
+ {
+ // If the DOF is constrained
+ DofConstraints::const_iterator
+ pos = _dof_constraints.find(elem_dofs[i]);
+
+ libmesh_assert (pos != _dof_constraints.end());
+
+ const DofConstraintRow& constraint_row = pos->second.first;
+
+// p refinement creates empty constraint rows
+// libmesh_assert (!constraint_row.empty());
+
+ for (DofConstraintRow::const_iterator
+ it=constraint_row.begin(); it != constraint_row.end();
+ ++it)
+ for (unsigned int j=0; j<new_elem_dofs.size(); j++)
+ if (new_elem_dofs[j] == it->first)
+ C(i,j) = it->second;
+ }
+ else
+ {
+ for (unsigned int j=0; j<new_elem_dofs.size(); j++)
+ if (new_elem_dofs[j] == elem_dofs[i])
+ C(i,j) = 1.;
+ }
+
+ // May need to do this recursively. It is possible
+ // that we just replaced a constrained DOF with another
+ // constrained DOF.
+ elem_dofs = new_elem_dofs;
+
+ DenseMatrix<Number> Cnew;
+ DenseVector<Number> Hnew;
+
+ this->build_constraint_matrix_and_vector (Cnew, Hnew, elem_dofs, true);
+
+ // If x = Cy + h and y = Dz + g
+ // Then x = (CD)z + (Cg + h)
+ C.vector_mult_add(H, 1, Hnew);
+
+ if ((C.n() == Cnew.m()) && // If the constraint matrix
+ (Cnew.n() == elem_dofs.size())) // is constrained...
+ C.right_multiply(Cnew);
+
+ libmesh_assert (C.n() == elem_dofs.size());
+ }
+
+ if (!called_recursively)
+ STOP_LOG("build_constraint_matrix_and_vector()", "DofMap");
+}
+
+
// NodeConstraints can fail in parallel when we try to look up a node
// that our processor doesn't have. Until we have a fix for that,
// let's just disable the allgather attempt.
@@ -1068,12 +1298,13 @@
libMesh::processor_id() - p) %
libMesh::n_processors();
- // Pack the dof constraint rows to push to procup
+ // Pack the dof constraint rows and rhs's to push to procup
std::vector<std::vector<unsigned int> > pushed_keys(pushed_ids[procup].size());
std::vector<std::vector<Real> > pushed_vals(pushed_ids[procup].size());
+ std::vector<Number> pushed_rhss(pushed_ids[procup].size());
for (unsigned int i = 0; i != pushed_ids[procup].size(); ++i)
{
- DofConstraintRow &row = _dof_constraints[pushed_ids[procup][i]];
+ DofConstraintRow &row = _dof_constraints[pushed_ids[procup][i]].first;
unsigned int row_size = row.size();
pushed_keys[i].reserve(row_size);
pushed_vals[i].reserve(row_size);
@@ -1083,16 +1314,18 @@
pushed_keys[i].push_back(j->first);
pushed_vals[i].push_back(j->second);
}
+ pushed_rhss[i] = _dof_constraints[pushed_ids[procup][i]].second;
}
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
// Pack the node constraint rows to push to procup
std::vector<std::vector<unsigned int> > pushed_node_keys(pushed_node_ids[procup].size());
std::vector<std::vector<Real> > pushed_node_vals(pushed_node_ids[procup].size());
+ std::vector<Point> pushed_node_offsets(pushed_node_ids[procup].size());
for (unsigned int i = 0; i != pushed_node_ids[procup].size(); ++i)
{
const Node *node = mesh.node_ptr(pushed_node_ids[procup][i]);
- NodeConstraintRow &row = _node_constraints[node];
+ NodeConstraintRow &row = _node_constraints[node].first;
unsigned int row_size = row.size();
pushed_node_keys[i].reserve(row_size);
pushed_node_vals[i].reserve(row_size);
@@ -1102,6 +1335,7 @@
pushed_node_keys[i].push_back(j->first->id());
pushed_node_vals[i].push_back(j->second);
}
+ pushed_node_offsets[i] = _node_constraints[node].second;
}
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
@@ -1109,28 +1343,36 @@
std::vector<unsigned int> pushed_ids_to_me;
std::vector<std::vector<unsigned int> > pushed_keys_to_me;
std::vector<std::vector<Real> > pushed_vals_to_me;
+ std::vector<Number> pushed_rhss_to_me;
Parallel::send_receive(procup, pushed_ids[procup],
procdown, pushed_ids_to_me);
Parallel::send_receive(procup, pushed_keys,
procdown, pushed_keys_to_me);
Parallel::send_receive(procup, pushed_vals,
procdown, pushed_vals_to_me);
+ Parallel::send_receive(procup, pushed_rhss,
+ procdown, pushed_rhss_to_me);
libmesh_assert (pushed_ids_to_me.size() == pushed_keys_to_me.size());
libmesh_assert (pushed_ids_to_me.size() == pushed_vals_to_me.size());
+ libmesh_assert (pushed_ids_to_me.size() == pushed_rhss_to_me.size());
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
// Trade pushed node constraint rows
std::vector<unsigned int> pushed_node_ids_to_me;
std::vector<std::vector<unsigned int> > pushed_node_keys_to_me;
std::vector<std::vector<Real> > pushed_node_vals_to_me;
+ std::vector<Point> pushed_node_offsets_to_me;
Parallel::send_receive(procup, pushed_node_ids[procup],
procdown, pushed_node_ids_to_me);
Parallel::send_receive(procup, pushed_node_keys,
procdown, pushed_node_keys_to_me);
Parallel::send_receive(procup, pushed_node_vals,
procdown, pushed_node_vals_to_me);
+ Parallel::send_receive(procup, pushed_node_offsets,
+ procdown, pushed_node_offsets_to_me);
libmesh_assert (pushed_node_ids_to_me.size() == pushed_node_keys_to_me.size());
libmesh_assert (pushed_node_ids_to_me.size() == pushed_node_vals_to_me.size());
+ libmesh_assert (pushed_node_ids_to_me.size() == pushed_node_offsets_to_me.size());
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
// Add the dof constraints that I've been sent
@@ -1144,11 +1386,12 @@
// add the one we were sent
if (!this->is_constrained_dof(constrained))
{
- DofConstraintRow &row = _dof_constraints[constrained];
+ DofConstraintRow &row = _dof_constraints[constrained].first;
for (unsigned int j = 0; j != pushed_keys_to_me[i].size(); ++j)
{
row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
}
+ _dof_constraints[constrained].second = pushed_rhss_to_me[i];
}
}
@@ -1165,13 +1408,14 @@
const Node *constrained = mesh.node_ptr(constrained_id);
if (!this->is_constrained_node(constrained))
{
- NodeConstraintRow &row = _node_constraints[constrained];
+ NodeConstraintRow &row = _node_constraints[constrained].first;
for (unsigned int j = 0; j != pushed_node_keys_to_me[i].size(); ++j)
{
const Node *key_node = mesh.node_ptr(pushed_node_keys_to_me[i][j]);
libmesh_assert(key_node);
row[key_node] = pushed_node_vals_to_me[i][j];
}
+ _node_constraints[constrained].second = pushed_node_offsets_to_me[i];
}
}
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
@@ -1228,7 +1472,7 @@
for (DoF_RCSet::iterator i = unexpanded_dofs.begin();
i != unexpanded_dofs.end(); ++i)
{
- DofConstraintRow &row = _dof_constraints[*i];
+ DofConstraintRow &row = _dof_constraints[*i].first;
for (DofConstraintRow::iterator j = row.begin();
j != row.end(); ++j)
dof_request_set.insert(j->first);
@@ -1238,7 +1482,7 @@
for (Node_RCSet::iterator i = unexpanded_nodes.begin();
i != unexpanded_nodes.end(); ++i)
{
- NodeConstraintRow &row = _node_constraints[*i];
+ NodeConstraintRow &row = _node_constraints[*i].first;
for (NodeConstraintRow::iterator j = row.begin();
j != row.end(); ++j)
{
@@ -1314,12 +1558,14 @@
node_row_keys(node_request_to_fill.size());
std::vector<std::vector<Real> > dof_row_vals(dof_request_to_fill.size()),
node_row_vals(node_request_to_fill.size());
+ std::vector<Number> dof_row_rhss(dof_request_to_fill.size());
+ std::vector<Point> node_row_rhss(node_request_to_fill.size());
for (unsigned int i=0; i != dof_request_to_fill.size(); ++i)
{
unsigned int constrained = dof_request_to_fill[i];
if (_dof_constraints.count(constrained))
{
- DofConstraintRow &row = _dof_constraints[constrained];
+ DofConstraintRow &row = _dof_constraints[constrained].first;
unsigned int row_size = row.size();
dof_row_keys[i].reserve(row_size);
dof_row_vals[i].reserve(row_size);
@@ -1329,6 +1575,7 @@
dof_row_keys[i].push_back(j->first);
dof_row_vals[i].push_back(j->second);
}
+ dof_row_rhss[i] = _dof_constraints[constrained].second;
}
}
@@ -1339,7 +1586,7 @@
const Node *constrained_node = mesh.node_ptr(constrained_id);
if (_node_constraints.count(constrained_node))
{
- const NodeConstraintRow &row = _node_constraints[constrained_node];
+ const NodeConstraintRow &row = _node_constraints[constrained_node].first;
unsigned int row_size = row.size();
node_row_keys[i].reserve(row_size);
node_row_vals[i].reserve(row_size);
@@ -1349,6 +1596,7 @@
node_row_keys[i].push_back(j->first->id());
node_row_vals[i].push_back(j->second);
}
+ node_row_rhss[i] = _node_constraints[constrained_node].second;
}
}
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
@@ -1358,6 +1606,8 @@
node_filled_keys;
std::vector<std::vector<Real> > dof_filled_vals,
node_filled_vals;
+ std::vector<Number> dof_filled_rhss;
+ std::vector<Point> node_filled_rhss;
Parallel::send_receive(procdown, dof_row_keys,
procup, dof_filled_keys);
Parallel::send_receive(procdown, dof_row_vals,
@@ -1367,11 +1617,15 @@
procup, node_filled_keys);
Parallel::send_receive(procdown, node_row_vals,
procup, node_filled_vals);
+ Parallel::send_receive(procdown, node_row_rhss,
+ procup, node_filled_rhss);
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
libmesh_assert (dof_filled_keys.size() == requested_dof_ids[procup].size());
libmesh_assert (dof_filled_vals.size() == requested_dof_ids[procup].size());
+ libmesh_assert (dof_filled_rhss.size() == requested_dof_ids[procup].size());
libmesh_assert (node_filled_keys.size() == requested_node_ids[procup].size());
libmesh_assert (node_filled_vals.size() == requested_node_ids[procup].size());
+ libmesh_assert (node_filled_rhss.size() == requested_node_ids[procup].size());
// Add any new constraint rows we've found
for (unsigned int i=0; i != requested_dof_ids[procup].size(); ++i)
@@ -1381,9 +1635,10 @@
if (!dof_filled_keys[i].empty())
{
unsigned int constrained = requested_dof_ids[procup][i];
- DofConstraintRow &row = _dof_constraints[constrained];
+ DofConstraintRow &row = _dof_constraints[constrained].first;
for (unsigned int j = 0; j != dof_filled_keys[i].size(); ++j)
row[dof_filled_keys[i][j]] = dof_filled_vals[i][j];
+ _dof_constraints[constrained].second = dof_filled_rhss[i];
// And prepare to check for more recursive constraints
unexpanded_dofs.insert(constrained);
@@ -1398,7 +1653,7 @@
{
unsigned int constrained_id = requested_node_ids[procup][i];
const Node* constrained_node = mesh.node_ptr(constrained_id);
- NodeConstraintRow &row = _node_constraints[constrained_node];
+ NodeConstraintRow &row = _node_constraints[constrained_node].first;
for (unsigned int j = 0; j != node_filled_keys[i].size(); ++j)
{
const Node* key_node =
@@ -1406,6 +1661,7 @@
libmesh_assert(key_node);
row[key_node] = node_filled_vals[i][j];
}
+ _node_constraints[constrained_node].second = node_filled_rhss[i];
// And prepare to check for more recursive constraints
unexpanded_nodes.insert(constrained_node);
@@ -1444,14 +1700,14 @@
libmesh_assert (pos != _dof_constraints.end());
- DofConstraintRow& constraint_row = pos->second;
+ DofConstraintRow& constraint_row = pos->second.first;
std::vector<unsigned int> constraints_to_expand;
for (DofConstraintRow::const_iterator
it=constraint_row.begin(); it != constraint_row.end();
++it)
- if (this->is_constrained_dof(it->first))
+ if (it->first != *i && this->is_constrained_dof(it->first))
{
unexpanded_set.insert(it->first);
constraints_to_expand.push_back(it->first);
@@ -1467,7 +1723,7 @@
libmesh_assert (subpos != _dof_constraints.end());
- const DofConstraintRow& subconstraint_row = subpos->second;
+ const DofConstraintRow& subconstraint_row = subpos->second.first;
for (DofConstraintRow::const_iterator
it=subconstraint_row.begin();
@@ -1518,7 +1774,7 @@
constrained_dof >= this->end_dof())
continue;
- DofConstraintRow& constraint_row = i->second;
+ DofConstraintRow& constraint_row = i->second.first;
for (DofConstraintRow::const_iterator
j=constraint_row.begin(); j != constraint_row.end();
++j)
@@ -1582,7 +1838,10 @@
// Add "this is zero" constraint rows for high p vertex
// dofs
for (unsigned int i = low_nc; i != high_nc; ++i)
- _dof_constraints[node->dof_number(sys_num,var,i)].clear();
+ {
+ _dof_constraints[node->dof_number(sys_num,var,i)].first.clear();
+ _dof_constraints[node->dof_number(sys_num,var,i)].second = 0.;
+ }
}
else
{
@@ -1593,7 +1852,8 @@
for (unsigned int j = low_nc; j != high_nc; ++j)
{
const unsigned int i = total_dofs - j - 1;
- _dof_constraints[node->dof_number(sys_num,var,i)].clear();
+ _dof_constraints[node->dof_number(sys_num,var,i)].first.clear();
+ _dof_constraints[node->dof_number(sys_num,var,i)].second = 0.;
}
}
}
------------------------------------------------------------------------------
Try before you buy = See our experts in action!
The most comprehensive online learning library for Microsoft developers
is just $99.99! Visual Studio, SharePoint, SQL - plus HTML5, CSS3, MVC3,
Metro Style Apps, more. Free future releases when you subscribe now!
http://p.sf.net/sfu/learndevnow-dev2
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
