On 4/28/10 1:38 PM, Roy Stogner wrote:
On Wed, 28 Apr 2010, Boyce Griffith wrote:
On 4/28/10 1:26 PM, Roy Stogner wrote:
On Wed, 28 Apr 2010, Boyce Griffith wrote:
For AMR meshes, should these sorts of constraints be applied only for
the nodes of level 0 elements? Or is this something that I don't have
to worry about explicitly?
You need to apply your constraints for all dofs that aren't already
constrained - everything except hanging nodes in adaptive h refinement
and their equivalent high-p dofs in adaptive p refinement. In your
case it shouldn't matter whether or not you constrain hanging nodes as
well, but if you do then turn off forbid_constraint_overwrite when
adding them; by default we assume adding two different constraints to
the same node is a bug.
Looking ahead to inhomogenous boundary conditions ---
A patch for that would be very appreciated, even if it only worked for
LAGRANGE variables at first, as long as the API was general enough to
be extensible later. This has *forever* been on our list of ideas
that are important/interesting enough to design but not quite
important/urgent enough to implement.
So, uhh, merely 7 months after bringing this up, here is a patch that
adds support for inhomogeneous DOF constraints. I think that most of it
is pretty straightforward. Basically, I changed the DofConstraints
typedef to keep track of both the rows of the constraint matrix, and a
right-hand-side value, which defaults to zero for the existing version
of add_constraint_row(), and I added a new version of
add_constraint_row() that allows for the specification of a RHS value.
Of course, these could be rolled into a single function.
The main gotcha, which took me a while to track down, is that
process_constraints() was not expecting a constrained DOF to be
"constrained to itself." I added a check in process_constraints() to
prevent such DOFs from being added to the list of expandable BCs ---
otherwise, the self-constraint winds up being removed from the row of
the constraint equation. Similar checks may be needed in
allgather_recursive_constraints(), but I have not added them, in part
because I don't know how to test this function (ParallelMesh doesn't
work yet, right?).
I've verified that ex4 and ex10 still work, both in serial and in
parallel. I've also made a version of ex4 that uses the DOF constraints
to impose Dirichlet BCs, which is attached, and which I've verified to
work in serial and in parallel.
If this seems useful, there is at least one issue that may warrant some
API changes:
The patched code does the "right thing" only when
asymmetric_constraint_rows=true when calling
DofMap::constrain_element_{matrix,vector,matrix_and_vector}(). This is
the default, but if asymmetric_constraint_rows=false, the code will
silently do the "wrong thing," probably yielding bizarre results. One
approach would be to keep track of whether a particular constraint row
is provided by the user or not, and to require
asymmetric_constraint_rows==true when such constraints are present.
-- Boyce
/* $Id: ex4.C 3894 2010-08-06 22:42:31Z roystgnr $ */
/* The Next Great Finite Element Library. */
/* Copyright (C) 2003 Benjamin S. Kirk */
/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */
/* License as published by the Free Software Foundation; either */
/* version 2.1 of the License, or (at your option) any later version. */
/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */
/* Lesser General Public License for more details. */
/* You should have received a copy of the GNU Lesser General Public */
/* License along with this library; if not, write to the Free Software */
/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
// <h1>Example 4 - Solving a 1D, 2D or 3D Poisson Problem in Parallel</h1>
//
// This is the fourth example program. It builds on
// the third example program by showing how to formulate
// the code in a dimension-independent way. Very minor
// changes to the example will allow the problem to be
// solved in one, two or three dimensions.
//
// This example will also introduce the PerfLog class
// as a way to monitor your code's performance. We will
// use it to instrument the matrix assembly code and look
// for bottlenecks where we should focus optimization efforts.
//
// This example also shows how to extend example 3 to run in
// parallel. Notice how litte has changed! The significant
// differences are marked with "PARALLEL CHANGE".
//
// NOTE: This version of ex4 imposes Dirichlet boundary
// conditions via inhomogeneous DOF constraints.
// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <math.h>
// Basic include file needed for the mesh functionality.
#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "exodusII_io.h"
#include "gnuplot_io.h"
#include "linear_implicit_system.h"
#include "equation_systems.h"
#include "boundary_info.h"
#include "exact_solution.h"
// Define the Finite Element object.
#include "fe.h"
// Define Gauss quadrature rules.
#include "quadrature_gauss.h"
// Define the DofMap, which handles degree of freedom
// indexing.
#include "dof_map.h"
// Define useful datatypes for finite element
// matrix and vector components.
#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "dense_matrix.h"
#include "dense_vector.h"
// Define the PerfLog, a performance logging utility.
// It is useful for timing events in a code and giving
// you an idea where bottlenecks lie.
#include "perf_log.h"
// The definition of a geometric element
#include "elem.h"
#include "string_to_enum.h"
#include "getpot.h"
// Bring in everything from the libMesh namespace
using namespace libMesh;
// Function prototype. This is the function that will assemble
// the linear system for our Poisson problem. Note that the
// function will take the \p EquationSystems object and the
// name of the system we are assembling as input. From the
// \p EquationSystems object we have acess to the \p Mesh and
// other objects we might need.
void assemble_poisson(EquationSystems& es,
const std::string& system_name);
// Function prototype. This is the function that will set DOF constraints to
// impose Dirichlet boundary conditions.
void constrain_poisson(EquationSystems& es,
const std::string& system_name);
// Exact solution function prototype.
Real exact_solution (const Real x,
const Real y = 0.,
const Real z = 0.);
// Version of exact_solution suitable to pass to an ExactSolution object, which
// is used to compute error norms.
Real exact_solution (const Point& p,
const Parameters&,
const std::string&,
const std::string&)
{
return exact_solution(p(0), p(1), p(2));
}
// Begin the main program.
int main (int argc, char** argv)
{
// Initialize libMesh and any dependent libaries, like in example 2.
LibMeshInit init (argc, argv);
// Declare a performance log for the main program
// PerfLog perf_main("Main Program");
// Create a GetPot object to parse the command line
GetPot command_line (argc, argv);
// Check for proper calling arguments.
if (argc < 3)
{
if (libMesh::processor_id() == 0)
std::cerr << "Usage:\n"
<<"\t " << argv[0] << " -d 2(3)" << " -n 15"
<< std::endl;
// This handy function will print the file name, line number,
// and then abort. Currrently the library does not use C++
// exception handling.
libmesh_error();
}
// Brief message to the user regarding the program name
// and command line arguments.
else
{
std::cout << "Running " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
}
// Read problem dimension from command line. Use int
// instead of unsigned since the GetPot overload is ambiguous
// otherwise.
int dim = 2;
if ( command_line.search(1, "-d") )
dim = command_line.next(dim);
// Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
// Create a mesh with user-defined dimension.
// Read number of elements from command line
int ps = 15;
if ( command_line.search(1, "-n") )
ps = command_line.next(ps);
// Read FE order from command line
std::string order = "SECOND";
if ( command_line.search(2, "-Order", "-o") )
order = command_line.next(order);
// Read FE Family from command line
std::string family = "LAGRANGE";
if ( command_line.search(2, "-FEFamily", "-f") )
family = command_line.next(family);
// Cannot use discontinuous basis.
if ((family == "MONOMIAL") || (family == "XYZ"))
{
if (libMesh::processor_id() == 0)
std::cerr << "ex4 currently requires a C^0 (or higher) FE basis." <<
std::endl;
libmesh_error();
}
Mesh mesh;
// Use the MeshTools::Generation mesh generator to create a uniform
// grid on the square [-1,1]^D. We instruct the mesh generator
// to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27
// elements in 3D. Building these higher-order elements allows
// us to use higher-order approximation, as in example 3.
Real halfwidth = dim > 1 ? 1. : 0.;
Real halfheight = dim > 2 ? 1. : 0.;
if ((family == "LAGRANGE") && (order == "FIRST"))
{
// No reason to use high-order geometric elements if we are
// solving with low-order finite elements.
MeshTools::Generation::build_cube (mesh,
ps,
(dim>1) ? ps : 0,
(dim>2) ? ps : 0,
-1., 1.,
-halfwidth, halfwidth,
-halfheight, halfheight,
(dim==1) ? EDGE2 :
((dim == 2) ? QUAD4 : HEX8));
}
else
{
MeshTools::Generation::build_cube (mesh,
ps,
(dim>1) ? ps : 0,
(dim>2) ? ps : 0,
-1., 1.,
-halfwidth, halfwidth,
-halfheight, halfheight,
(dim==1) ? EDGE3 :
((dim == 2) ? QUAD9 : HEX27));
}
// Print information about the mesh to the screen.
mesh.print_info();
// Create an equation systems object.
EquationSystems equation_systems (mesh);
// Declare the system and its variables.
// Create a system named "Poisson"
LinearImplicitSystem& system =
equation_systems.add_system<LinearImplicitSystem> ("Poisson");
// Add the variable "u" to "Poisson". "u"
// will be approximated using second-order approximation.
system.add_variable("u",
Utility::string_to_enum<Order> (order),
Utility::string_to_enum<FEFamily>(family));
// Give the system a pointer to the matrix assembly
// function.
system.attach_assemble_function (assemble_poisson);
// Give the system a pointer to the DOF constraint function.
system.attach_constraint_function (constrain_poisson);
// Initialize the data structures for the equation system.
equation_systems.init();
// Print information about the system to the screen.
equation_systems.print_info();
mesh.print_info();
// Solve the system "Poisson", just like example 2.
equation_systems.get_system("Poisson").solve();
// Evaluate the error in the solution.
ExactSolution exact_soln(equation_systems);
exact_soln.attach_exact_value(exact_solution);
exact_soln.compute_error("Poisson", "u");
const double l2_err = exact_soln.l2_error("Poisson", "u");
const double l1_err = exact_soln.l1_error("Poisson", "u");
const double l_inf_err = exact_soln.l_inf_error("Poisson", "u");
if (libMesh::processor_id() == 0)
std::cout << "l2 error: " << l2_err << "\n"
<< "l1 error: " << l1_err << "\n"
<< "l_inf error: " << l_inf_err << "\n";
const DofMap& dof_map = system.get_dof_map();
std::pair<Real,Real> err = dof_map.max_constraint_error(system);
if (libMesh::processor_id() == 0)
std::cout << "max constraint abs error: " << err.first << "\n"
<< "max constraint rel error: " << err.second << "\n";
// After solving the system write the solution
// to a GMV-formatted plot file.
if(dim == 1)
{
GnuPlotIO plot(mesh,"Example 4, 1D",GnuPlotIO::GRID_ON);
plot.write_equation_systems("out_1",equation_systems);
}
#ifdef LIBMESH_HAVE_EXODUS_API
else
{
ExodusII_IO (mesh).write_equation_systems ((dim == 3) ?
"out_3.exd" : "out_2.exd",equation_systems);
}
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
// All done.
return 0;
}
//
//
//
// We now define the matrix assembly function for the Poisson system. We need
// to compute element matrices and right-hand sides. Inhomogeneous Dirichlet
// boundary conditions are handled by the DofMap.
void assemble_poisson(EquationSystems& es,
const std::string& system_name)
{
// It is a good idea to make sure we are assembling
// the proper system.
libmesh_assert (system_name == "Poisson");
// Declare a performance log. Give it a descriptive
// string to identify what part of the code we are
// logging, since there may be many PerfLogs in an
// application.
PerfLog perf_log ("Matrix Assembly");
// Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
// The dimension that we are running
const unsigned int dim = mesh.mesh_dimension();
// Get a reference to the LinearImplicitSystem we are solving
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Poisson");
// A reference to the \p DofMap object for this system. The \p DofMap
// object handles the index translation from node and element numbers
// to degree of freedom numbers. We will talk more about the \p DofMap
// in future examples.
const DofMap& dof_map = system.get_dof_map();
// Get a constant reference to the Finite Element type
// for the first (and only) variable in the system.
FEType fe_type = dof_map.variable_type(0);
// Build a Finite Element object of the specified type. Since the
// \p FEBase::build() member dynamically creates memory we will
// store the object as an \p AutoPtr<FEBase>. This can be thought
// of as a pointer that will clean up after itself.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
// A 5th order Gauss quadrature rule for numerical integration.
QGauss qrule (dim, FIFTH);
// Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
// Declare a special finite element object for
// boundary integration.
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
// Boundary integration requires one quadraure rule,
// with dimensionality one less than the dimensionality
// of the element.
QGauss qface(dim-1, FIFTH);
// Tell the finte element object to use our
// quadrature rule.
fe_face->attach_quadrature_rule (&qface);
// Here we define some references to cell-specific data that
// will be used to assemble the linear system.
// We begin with the element Jacobian * quadrature weight at each
// integration point.
const std::vector<Real>& JxW = fe->get_JxW();
// The physical XY locations of the quadrature points on the element.
// These might be useful for evaluating spatially varying material
// properties at the quadrature points.
const std::vector<Point>& q_point = fe->get_xyz();
// The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
// The element shape function gradients evaluated at the quadrature
// points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
// Define data structures to contain the element matrix
// and right-hand-side vector contribution. Following
// basic finite element terminology we will denote these
// "Ke" and "Fe". More detail is in example 3.
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
// This vector will hold the degree of freedom indices for
// the element. These define where in the global system
// the element degrees of freedom get mapped.
std::vector<unsigned int> dof_indices;
// Now we will loop over all the elements in the mesh.
// We will compute the element matrix and right-hand-side
// contribution. See example 3 for a discussion of the
// element iterators.
MeshBase::const_element_iterator el =
mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
// Start logging the shape function initialization.
// This is done through a simple function call with
// the name of the event to log.
perf_log.push("elem init");
// Store a pointer to the element we are currently
// working on. This allows for nicer syntax later.
const Elem* elem = *el;
// Get the degree of freedom indices for the
// current element. These define where in the global
// matrix and right-hand-side this element will
// contribute to.
dof_map.dof_indices (elem, dof_indices);
// Compute the element-specific data for the current
// element. This involves computing the location of the
// quadrature points (q_point) and the shape functions
// (phi, dphi) for the current element.
fe->reinit (elem);
// Zero the element matrix and right-hand side before
// summing them. We use the resize member here because
// the number of degrees of freedom might have changed from
// the last element. Note that this will be the case if the
// element type is different (i.e. the last element was a
// triangle, now we are on a quadrilateral).
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
// Stop logging the shape function initialization.
// If you forget to stop logging an event the PerfLog
// object will probably catch the error and abort.
perf_log.pop("elem init");
// Now we will build the element matrix. This involves
// a double loop to integrate the test funcions (i) against
// the trial functions (j).
//
// We have split the numeric integration into two loops
// so that we can log the matrix and right-hand-side
// computation seperately.
//
// Now start logging the element matrix computation
perf_log.push ("Ke");
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
// Stop logging the matrix computation
perf_log.pop ("Ke");
// Now we build the element right-hand-side contribution.
// This involves a single loop in which we integrate the
// "forcing function" in the PDE against the test functions.
//
// Start logging the right-hand-side computation
perf_log.push ("Fe");
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
// fxy is the forcing function for the Poisson equation.
// In this case we set fxy to be a finite difference
// Laplacian approximation to the (known) exact solution.
//
// We will use the second-order accurate FD Laplacian
// approximation, which in 2D on a structured grid is
//
// u_xx + u_yy = (u(i-1,j) + u(i+1,j) +
// u(i,j-1) + u(i,j+1) +
// -4*u(i,j))/h^2
//
// Since the value of the forcing function depends only
// on the location of the quadrature point (q_point[qp])
// we will compute it here, outside of the i-loop
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
const Real z = q_point[qp](2);
const Real eps = 1.e-3;
const Real uxx = (exact_solution(x-eps,y,z) +
exact_solution(x+eps,y,z) +
-2.*exact_solution(x,y,z))/eps/eps;
const Real uyy = (exact_solution(x,y-eps,z) +
exact_solution(x,y+eps,z) +
-2.*exact_solution(x,y,z))/eps/eps;
const Real uzz = (exact_solution(x,y,z-eps) +
exact_solution(x,y,z+eps) +
-2.*exact_solution(x,y,z))/eps/eps;
Real fxy;
if(dim==1)
{
// In 1D, compute the rhs by differentiating the
// exact solution twice.
const Real pi = libMesh::pi;
fxy = (0.25*pi*pi)*sin(.5*pi*x);
}
else
{
fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));
}
// Add the RHS contribution
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*fxy*phi[i][qp];
}
// Stop logging the right-hand-side computation
perf_log.pop ("Fe");
// If this assembly program were to be used on an adaptive mesh,
// we would have to apply any hanging node constraint equations
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
// The element matrix and right-hand-side are now built
// for this element. Add them to the global matrix and
// right-hand-side vector. The \p SparseMatrix::add_matrix()
// and \p NumericVector::add_vector() members do this for us.
// Start logging the insertion of the local (element)
// matrix and vector into the global matrix and vector
perf_log.push ("matrix insertion");
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
// Start logging the insertion of the local (element)
// matrix and vector into the global matrix and vector
perf_log.pop ("matrix insertion");
}
// That's it. We don't need to do anything else to the
// PerfLog. When it goes out of scope (at this function return)
// it will print its log to the screen. Pretty easy, huh?
}
// Function to set inhomogeneous DOF constraints, which are used to impose
// Dirichlet boundary conditions.
void constrain_poisson(EquationSystems& es,
const std::string& system_name)
{
// It is a good idea to make sure we are assembling
// the proper system.
libmesh_assert (system_name == "Poisson");
// Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
// Get a reference to the LinearImplicitSystem we are solving.
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Poisson");
const unsigned int system_number = system.number();
const unsigned int var_number = system.variable_number("u");
// Setup Dirichlet boundary conditions as DOF constraints. (This is not
// especially elegant, but does the job.)
DofMap& dof_map = system.get_dof_map();
MeshBase::const_element_iterator el =
mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
const Elem* const elem = *el;
for (unsigned int side=0; side<elem->n_sides(); side++)
if (elem->neighbor(side) == NULL)
{
for (unsigned node=0; node<elem->n_nodes(); ++node)
{
if (elem->is_node_on_side(node,side))
{
const Node& n = *(elem->get_node(node));
const unsigned int dof_id = n.dof_number(system_number,
var_number, 0);
const Real x = n(0);
const Real y = n(1);
const Real z = n(2);
DofConstraintRow row;
row.insert(std::make_pair(dof_id,1.0));
const Real rhs = exact_solution(x,y,z);
dof_map.add_constraint_row(dof_id, row, rhs, false);
}
}
}
}
}
Index: include/base/dof_map.h
===================================================================
--- include/base/dof_map.h (revision 4138)
+++ include/base/dof_map.h (working copy)
@@ -80,7 +80,7 @@
* 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> > >
{
@@ -355,12 +355,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 constraint row
@@ -373,7 +384,7 @@
*/
DofConstraints::const_iterator constraint_rows_end() const
{ return _dof_constraints.end(); }
-
+
/**
* @returns true if the degree of freedom dof is constrained,
* false otherwise.
Index: src/base/dof_map.C
===================================================================
--- src/base/dof_map.C (revision 4138)
+++ src/base/dof_map.C (working copy)
@@ -1854,7 +1854,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]
Index: src/base/dof_map_constraints.C
===================================================================
--- src/base/dof_map_constraints.C (revision 4138)
+++ src/base/dof_map_constraints.C (working copy)
@@ -191,6 +191,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.
@@ -202,7 +203,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);
}
@@ -218,7 +219,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";
@@ -228,6 +230,8 @@
os << " (" << pos->first << ","
<< pos->second << ")\t";
+ os << "rhs: " << rhs;
+
os << std::endl;
}
}
@@ -283,7 +287,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());
@@ -356,7 +360,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());
@@ -390,7 +394,12 @@
if (this->is_constrained_dof(elem_dofs[i]))
{
// If the DOF is constrained
- rhs(i) = 0.;
+ DofConstraints::const_iterator
+ pos = _dof_constraints.find(elem_dofs[i]);
+
+ libmesh_assert (pos != _dof_constraints.end());
+
+ rhs(i) = -pos->second.second;
}
} // end if is constrained...
@@ -461,7 +470,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());
@@ -519,7 +528,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.
@@ -576,7 +590,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.
@@ -671,9 +690,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)
@@ -747,7 +766,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)
@@ -797,7 +821,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());
@@ -849,7 +873,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());
@@ -946,12 +970,13 @@
libMesh::processor_id() - p) %
libMesh::n_processors();
- // Pack the constraint rows to push to procup
+ // Pack the 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);
@@ -961,20 +986,25 @@
pushed_keys[i].push_back(j->first);
pushed_vals[i].push_back(j->second);
}
+ pushed_rhss[i] = _dof_constraints[pushed_ids[procup][i]].second;
}
// Trade pushed constraint rows
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());
// Add the constraints that I've been sent
for (unsigned int i = 0; i != pushed_ids_to_me.size(); ++i)
@@ -987,11 +1017,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];
}
}
}
@@ -1031,7 +1062,7 @@
for (RCSet::iterator i = unexpanded_set.begin();
i != unexpanded_set.end(); ++i)
{
- DofConstraintRow &row = _dof_constraints[*i];
+ DofConstraintRow &row = _dof_constraints[*i].first;
for (DofConstraintRow::iterator j = row.begin();
j != row.end(); ++j)
request_set.insert(j->first);
@@ -1078,12 +1109,13 @@
// Fill those requests
std::vector<std::vector<unsigned int> >
row_keys(request_to_fill.size());
std::vector<std::vector<Real> > row_vals(request_to_fill.size());
+ std::vector<Number> row_rhss(request_to_fill.size());
for (unsigned int i=0; i != request_to_fill.size(); ++i)
{
unsigned int constrained = 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();
row_keys[i].reserve(row_size);
row_vals[i].reserve(row_size);
@@ -1093,18 +1125,23 @@
row_keys[i].push_back(j->first);
row_vals[i].push_back(j->second);
}
+ row_rhss[i] = _dof_constraints[constrained].second;
}
}
// Trade back the results
std::vector<std::vector<unsigned int> > filled_keys;
std::vector<std::vector<Real> > filled_vals;
+ std::vector<Number> filled_rhss;
Parallel::send_receive(procdown, row_keys,
procup, filled_keys);
Parallel::send_receive(procdown, row_vals,
procup, filled_vals);
+ Parallel::send_receive(procdown, row_rhss,
+ procup, filled_rhss);
libmesh_assert (filled_keys.size() == requested_ids[procup].size());
libmesh_assert (filled_vals.size() == requested_ids[procup].size());
+ libmesh_assert (filled_rhss.size() == requested_ids[procup].size());
// Add any new constraint rows we've found
for (unsigned int i=0; i != requested_ids[procup].size(); ++i)
@@ -1113,9 +1150,10 @@
if (!filled_keys[i].empty())
{
unsigned int constrained = requested_ids[procup][i];
- DofConstraintRow &row = _dof_constraints[constrained];
+ DofConstraintRow &row = _dof_constraints[constrained].first;
for (unsigned int j = 0; j != filled_keys[i].size(); ++j)
row[filled_keys[i][j]] = filled_vals[i][j];
+ _dof_constraints[constrained].second = filled_rhss[i];
// And prepare to check for more recursive constraints
unexpanded_set.insert(constrained);
@@ -1152,14 +1190,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);
@@ -1175,7 +1213,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();
@@ -1226,7 +1264,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)
@@ -1290,7 +1328,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
{
@@ -1301,7 +1342,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.;
}
}
}
Index: src/fe/fe_base.C
===================================================================
--- src/fe/fe_base.C (revision 4138)
+++ src/fe/fe_base.C (working copy)
@@ -1913,7 +1913,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));
@@ -2154,7 +2154,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)
{
@@ -2202,7 +2202,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));
Index: src/fe/fe_lagrange.C
===================================================================
--- src/fe/fe_lagrange.C (revision 4138)
+++ src/fe/fe_lagrange.C (working copy)
@@ -763,7 +763,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));
------------------------------------------------------------------------------
Increase Visibility of Your 3D Game App & Earn a Chance To Win $500!
Tap into the largest installed PC base & get more eyes on your game by
optimizing for Intel(R) Graphics Technology. Get started today with the
Intel(R) Software Partner Program. Five $500 cash prizes are up for grabs.
http://p.sf.net/sfu/intelisp-dev2dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
