On Thu, Oct 6, 2016 at 10:44 PM, David Knezevic <david.kneze...@akselos.com>
wrote:
> On Thu, Oct 6, 2016 at 10:34 PM, Roy Stogner <royst...@ices.utexas.edu>
> wrote:
>
>>
>>
>> On Thu, 6 Oct 2016, David Knezevic wrote:
>>
>> I'm using GhostingFunctor for a contact solve, in which I consider a 1/4
>>> domain with partial Dirichlet boundary conditions that impose a symmetry
>>> condition (i.e. displacement normal to the symmetry boundary is clamped
>>> to zero, and tangential displacement is unconstrained). This means that I
>>> have Dirichlet constraints that affect the dofs on the contact surface.
>>> What I find is that the solve works fine in serial, but in parallel the
>>> nonlinear convergence fails, presumably because of an incorrect Jacobian.
>>> I have actually run into this exact issue before (when I was augmenting
>>> the sparsity pattern "manually", prior to GhostingFunctor) and I found
>>> that the issue was that the dof constraints on the contact surface were
>>> not being communicated in parallel, which caused the incorrect Jacobian
>>> in parallel. I fixed it by adding artificial Edge2 elements into the
>>> mesh (like in systems_of_equations_ex8) to ensure that the dof constraints
>>> are communicated properly in parallel.
>>>
>>> So, my question is, is there a way to achieve the necessary dof
>>> constraint communication with the new GhostingFunctor API? I had hoped that
>>> using
>>> "add_algebraic_ghosting_functor" would do the job, but it apparently
>>> doesn't.
>>>
>>
>> Hmm... If you only needed algebraic ghosting, then
>> add_algebraic_ghosting_functor should have been sufficient -
>> processors won't know about all their ghosted dofs' constraints, but
>> the ghosted dofs will get properly added to the send_list, and
>> enforce_constraints_exactly() will ensure that the dofs, once
>> constrained on the processors which own them, will have their
>> values updated on all the processors which ghost them.
>>
>> But you need to augment the sparsity pattern too, so you should be
>> using add_coupling_functor instead... and in *that* case, we're
>> broken, aren't we? You build a Jacobian connecting the remotely
>> coupled dofs, and you try to hit it with constrain_element_foo() or
>> heterogeneously_constrain_element_foo(), but the processor isn't aware
>> of all the ghosted dofs' constraints, so the element constraint matrix
>> is wrong and so is your final answer.
>>
>
>
> Yep, that's exactly my understanding of the issue.
>
>
>
>> I think the proper fix is to call the coupling functors in
>> scatter_constraints(), then query the processors who own the elements
>> which are to be coupled for any constraints on those elements' dofs.
>> I can take a crack at that tomorrow or Monday. Any chance you could
>> set up a unit test (or even just stuff an assertion into the misc ex9
>> example?) that checks for the problem?
>>
>
> That'd be great, thanks! I'll be happy to try it out once it's ready. I'll
> also look into making a test case that can be added to libMesh.
>
I've attached a modified version of misc_ex9 that attaches constraints on
one side of the "crack" and checks if the dof constraints are present
during assembly. This passes in serial but fails in parallel because
constraints are not communicated on the GhostingFunctor-coupled-dofs.
The extra constraints I added mean that the problem doesn't make physical
sense anymore unfortunately, but at least it tests for the constraint
issue. Roy: I'm not sure if this is appropriate for a unit test, but it
should at least be helpful for when you want to check your implementation.
David
// The libMesh Finite Element Library.
// Copyright (C) 2002-2016 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
// 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>Miscellaneous Example 9 - Implement an interface term to model
// a thermal "film resistance"</h1>
// \author David Knezevic
// \date 2013
//
// In this example we solve a Poisson problem, -\Laplacian u = f, with
// a non-standard interface condition on the domain interior which
// models a thermal "film resistance". The interface condition
// requires continuity of flux, and a jump in temperature proportional
// to the flux:
// \nabla u_1 \cdot n = \nabla u_2 \cdot n,
// u_1 - u_2 = R * \nabla u \cdot n
//
// To implement this PDE, we use two mesh subdomains, \Omega_1 and
// \Omega_2, with coincident boundaries, but which are not connected
// in the FE sense. Let \Gamma denote the coincident boundary. The
// term on \Gamma takes the form:
//
// 1/R * \int_\Gamma (u_1 - u_2) (v_1 - v_2) ds,
//
// where u_1, u_2 (resp. v_1, v_2) are the trial (resp. test)
// functions on either side of \Gamma. We implement this condition
// using C0 basis functions, but the "crack" in the mesh at \Gamma
// permits a discontinuity in the solution. We also impose a heat flux
// on the bottom surface of the mesh, and a zero Dirichlet condition
// on the top surface.
//
// In order to implement the interface condition, we need to augment
// the matrix sparsity pattern, which is handled by the class
// AugmentSparsityPatternOnInterface.
// C++ include files that we need
#include <iostream>
#include <limits>
// libMesh includes
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/getpot.h"
#include "libmesh/elem.h"
#include "libmesh/fe_interface.h"
#include "libmesh/boundary_info.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/zero_function.h"
#include "libmesh/dirichlet_boundaries.h"
// local includes
#include "augment_sparsity_on_interface.h"
// define the boundary IDs in the mesh
#define MIN_Z_BOUNDARY 15
#define MAX_Z_BOUNDARY 0
#define CRACK_BOUNDARY_LOWER 10
#define CRACK_BOUNDARY_UPPER 5
// Bring in everything from the libMesh namespace
using namespace libMesh;
/**
* Assemble the system matrix and rhs vector.
*/
void assemble_poisson(EquationSystems & es,
const ElementSideMap & lower_to_upper);
// The main program.
int main (int argc, char ** argv)
{
// Initialize libMesh.
LibMeshInit init (argc, argv);
// This example uses an ExodusII input file
#ifndef LIBMESH_HAVE_EXODUS_API
libmesh_example_requires(false, "--enable-exodus");
#endif
// The sparsity augmentation code requires PETSc
libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
// Skip this 3D example if libMesh was compiled as 1D or 2D-only.
libmesh_example_requires(3 <= LIBMESH_DIM, "3D support");
GetPot command_line (argc, argv);
Real R = 2.;
if (command_line.search(1, "-R"))
R = command_line.next(R);
Mesh mesh(init.comm());
EquationSystems equation_systems (mesh);
LinearImplicitSystem & system =
equation_systems.add_system<LinearImplicitSystem> ("Poisson");
system.add_variable("u", FIRST, LAGRANGE);
// We want to call assemble_poisson "manually" so that we can pass in
// lower_to_upper, hence set assemble_before_solve = false
system.assemble_before_solve = false;
// Impose zero Dirichlet boundary condition
std::set<boundary_id_type> boundary_ids;
boundary_ids.insert(MAX_Z_BOUNDARY);
boundary_ids.insert(CRACK_BOUNDARY_UPPER);
std::vector<unsigned int> variables;
variables.push_back(0);
ZeroFunction<> zf;
DirichletBoundary dirichlet_bc(boundary_ids,
variables,
&zf);
system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
// Attach an object to the DofMap that will augment the sparsity pattern
// due to the degrees-of-freedom on the "crack"
//
// By attaching this object *before* reading our mesh, we also
// ensure that the connected elements will not be deleted on a
// distributed mesh.
AugmentSparsityOnInterface augment_sparsity
(mesh, CRACK_BOUNDARY_LOWER, CRACK_BOUNDARY_UPPER);
system.get_dof_map().add_coupling_functor(augment_sparsity);
// This example code has not been written to cope with a distributed mesh
{
Mesh mesh_1 (init.comm(), 3);
MeshTools::Generation::build_cube (mesh_1,
4,
4,
20,
0., 4.,
0., 4.,
0., 20.,
HEX8);
Mesh mesh_2 (init.comm(), 3);
MeshTools::Generation::build_cube (mesh_2,
4,
4,
20,
0., 4.,
0., 4.,
20., 40.,
HEX8);
for(boundary_id_type bc_id = 0; bc_id < 6; bc_id++)
{
// Convert the BC IDs in one of the compnents to prevent "overlap" of
// the boundary IDs.
MeshTools::Modification::change_boundary_id(mesh_2,
bc_id,
bc_id+10);
}
// Call stitch_meshes to merge the two meshes, but pass invalid_id
// so that we do not do any actuall stitching.
mesh_2.stitch_meshes(mesh_1, BoundaryInfo::invalid_id, BoundaryInfo::invalid_id);
mesh_2.write("two_beams.exo");
}
mesh.read("two_beams.exo");
equation_systems.init();
equation_systems.print_info();
// Set the jump term coefficient, it will be used in assemble_poisson
equation_systems.parameters.set<Real>("R") = R;
// Assemble and then solve
assemble_poisson(equation_systems,
augment_sparsity.get_lower_to_upper());
system.solve();
#ifdef LIBMESH_HAVE_EXODUS_API
// Plot the solution
ExodusII_IO (mesh).write_equation_systems ("solution.exo",
equation_systems);
#endif
#ifdef LIBMESH_ENABLE_AMR
// Possibly solve on a refined mesh next.
MeshRefinement mesh_refinement (mesh);
unsigned int n_refinements = 0;
if (command_line.search("-n_refinements"))
n_refinements = command_line.next(0);
for (unsigned int r = 0; r != n_refinements; ++r)
{
std::cout << "Refining the mesh" << std::endl;
mesh_refinement.uniformly_refine ();
equation_systems.reinit();
assemble_poisson(equation_systems,
augment_sparsity.get_lower_to_upper());
system.solve();
#ifdef LIBMESH_HAVE_EXODUS_API
// Plot the refined solution
std::ostringstream out;
out << "solution_" << r << ".exo";
ExodusII_IO (mesh).write_equation_systems (out.str(),
equation_systems);
#endif
}
#endif
return 0;
}
void assemble_poisson(EquationSystems & es,
const ElementSideMap & lower_to_upper)
{
const MeshBase & mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
Real R = es.parameters.get<Real>("R");
LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Poisson");
const DofMap & dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
UniquePtr<FEBase> fe_elem_face (FEBase::build(dim, fe_type));
UniquePtr<FEBase> fe_neighbor_face (FEBase::build(dim, fe_type));
QGauss qrule (dim, fe_type.default_quadrature_order());
QGauss qface(dim-1, fe_type.default_quadrature_order());
fe->attach_quadrature_rule (&qrule);
fe_elem_face->attach_quadrature_rule (&qface);
fe_neighbor_face->attach_quadrature_rule (&qface);
const std::vector<Real> & JxW = fe->get_JxW();
const std::vector<std::vector<Real> > & phi = fe->get_phi();
const std::vector<std::vector<RealGradient> > & dphi = fe->get_dphi();
const std::vector<Real> & JxW_face = fe_elem_face->get_JxW();
const std::vector<Point> & qface_points = fe_elem_face->get_xyz();
const std::vector<std::vector<Real> > & phi_face = fe_elem_face->get_phi();
const std::vector<std::vector<Real> > & phi_neighbor_face = fe_neighbor_face->get_phi();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
DenseMatrix<Number> Kne;
DenseMatrix<Number> Ken;
DenseMatrix<Number> Kee;
DenseMatrix<Number> Knn;
std::vector<dof_id_type> dof_indices;
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 * elem = *el;
dof_map.dof_indices (elem, dof_indices);
const unsigned int n_dofs = dof_indices.size();
fe->reinit (elem);
Ke.resize (n_dofs, n_dofs);
Fe.resize (n_dofs);
// Assemble element interior terms for the matrix
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
for (unsigned int i=0; i<n_dofs; i++)
for (unsigned int j=0; j<n_dofs; j++)
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
// Boundary flux provides forcing in this example
{
for (unsigned int side=0; side<elem->n_sides(); side++)
if (elem->neighbor_ptr(side) == libmesh_nullptr)
{
if (mesh.get_boundary_info().has_boundary_id (elem, side, MIN_Z_BOUNDARY))
{
fe_elem_face->reinit(elem, side);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW_face[qp] * phi_face[i][qp];
}
}
}
// Add boundary terms on the crack
{
for (unsigned int side=0; side<elem->n_sides(); side++)
if (elem->neighbor_ptr(side) == libmesh_nullptr)
{
// Found the lower side of the crack. Assemble terms due to lower and upper in here.
if (mesh.get_boundary_info().has_boundary_id (elem, side, CRACK_BOUNDARY_LOWER))
{
fe_elem_face->reinit(elem, side);
ElementSideMap::const_iterator ltu_it =
lower_to_upper.find(std::make_pair(elem, side));
libmesh_assert(ltu_it != lower_to_upper.end());
const Elem * neighbor = ltu_it->second;
std::vector<Point> qface_neighbor_points;
FEInterface::inverse_map (elem->dim(), fe->get_fe_type(),
neighbor, qface_points, qface_neighbor_points);
fe_neighbor_face->reinit(neighbor, &qface_neighbor_points);
std::vector<dof_id_type> neighbor_dof_indices;
dof_map.dof_indices (neighbor, neighbor_dof_indices);
const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
Kne.resize (n_neighbor_dofs, n_dofs);
Ken.resize (n_dofs, n_neighbor_dofs);
Kee.resize (n_dofs, n_dofs);
Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
// Lower-to-lower coupling term
for (unsigned int qp=0; qp<qface.n_points(); qp++)
for (unsigned int i=0; i<n_dofs; i++)
for (unsigned int j=0; j<n_dofs; j++)
Kee(i,j) -= JxW_face[qp] * (1./R)*(phi_face[i][qp] * phi_face[j][qp]);
// Lower-to-upper coupling term
for (unsigned int qp=0; qp<qface.n_points(); qp++)
for (unsigned int i=0; i<n_dofs; i++)
for (unsigned int j=0; j<n_neighbor_dofs; j++)
Ken(i,j) += JxW_face[qp] * (1./R)*(phi_face[i][qp] * phi_neighbor_face[j][qp]);
// Upper-to-upper coupling term
for (unsigned int qp=0; qp<qface.n_points(); qp++)
for (unsigned int i=0; i<n_neighbor_dofs; i++)
for (unsigned int j=0; j<n_neighbor_dofs; j++)
Knn(i,j) -= JxW_face[qp] * (1./R)*(phi_neighbor_face[i][qp] * phi_neighbor_face[j][qp]);
// Upper-to-lower coupling term
for (unsigned int qp=0; qp<qface.n_points(); qp++)
for (unsigned int i=0; i<n_neighbor_dofs; i++)
for (unsigned int j=0; j<n_dofs; j++)
Kne(i,j) += JxW_face[qp] * (1./R)*(phi_neighbor_face[i][qp] * phi_face[j][qp]);
system.matrix->add_matrix(Kne, neighbor_dof_indices, dof_indices);
system.matrix->add_matrix(Ken, dof_indices, neighbor_dof_indices);
system.matrix->add_matrix(Kee, dof_indices);
system.matrix->add_matrix(Knn, neighbor_dof_indices);
// Find the dofs on neighbor that are on CRACK_BOUNDARY_UPPER
for(unsigned int neighbor_side=0;
neighbor_side<neighbor->n_sides();
neighbor_side++)
{
if(mesh.get_boundary_info().has_boundary_id (
neighbor, neighbor_side, CRACK_BOUNDARY_UPPER))
{
for(unsigned int neighbor_node_index=0;
neighbor_node_index<neighbor->n_nodes();
neighbor_node_index++)
{
if(neighbor->is_node_on_side(neighbor_node_index, neighbor_side))
{
dof_id_type neighbor_node_id = neighbor->node(neighbor_node_index);
if(!dof_map.is_constrained_dof(neighbor_node_id))
{
libmesh_error_msg("Neighbor dofs should have constraints on them");
}
}
}
}
}
}
}
}
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
}
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, SlashDot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
