[Libmesh-users] Periodic boundary condition does not get correct result (with code)

[Shengli X. <she...@gm...>]

Micheal, Wout, Ben
I am testing a code for the problem of 2D plane stress. I want to
impose periodic boundary condition just like u1=3Du2. I mainly consider
the post of Ben and Micheal.
I don't konw where to call the subroutine for periodic bc. The result
considering periodic bc is the same as not considering periodic bc. I
don't konw how to do. Maybe The part of code considering periodic bc
is wrong.

I have mailed to Wout for this problem. Maybe I need to add some code to
system.h and system.C.

The source code is attached.
The example is: 1x1 domain, 10x10 mesh, Quad4, fixed bc on x=3D0, p=3D1 on
x=3D1, consider periodic bc of ux with node (0.5,0) and (0.8,0).

// 2D elastic code

/* 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  US=
A
*/

#include <iostream>
#include <fstream>
#include <algorithm>
#include <math.h>
#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "gmv_io.h"
#include "linear_implicit_system.h"
#include "equation_systems.h"
#include "fe.h"
#include "quadrature_gauss.h"

#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "dense_matrix.h"
#include "dense_vector.h"
// add by shengli
#include "elem.h"

#include "dof_map.h"
#include "perf_log.h"
// for assembling the element matrix and vector on a component-by-component
basis
#include "dense_submatrix.h"
//#include "dense_sub
void assemble_elastic(EquationSystems& es, const std::string& system_name);

// impose periodic boundary condition
void apply_periodic_bc(EquationSystems &es, const std::string &system_name)=
;

int main (int argc, char** argv)
{
  libMesh::init (argc, argv);
  {
    std::cout << "Running " << argv[0];
    for (int i=3D1; i<argc; i++)
      std::cout << " " << argv[i];
    std::cout << std::endl << std::endl;

    PerfLog perf("Main Program");
    perf.start_event("program init");
    Mesh mesh (2);
    // Use the MeshTools::Generation mesh generator to create a uniform
    // grid on the square [0,1]^2.
    MeshTools::Generation::build_square(mesh,10,10,0.,1.,0.,1.,QUAD4);

    //  MeshData mesh_data(mesh);
    //mesh_data.activate();
    //  mesh_data.enable_compatibility_mode();

    //  mesh.read("in.xda",&mesh_data);
    //  mesh_data.read("data.xta");
    mesh.print_info();

    //  EquationSystems equation_systems (mesh,&mesh_data);
    EquationSystems equation_systems (mesh);
    LinearImplicitSystem &system =3D
      equation_systems.add_system<LinearImplicitSystem>("elastic");

    // equation_systems.add_system<LinearImplicitSystem> ("elastic");
    //equation_systems("elastic").add_variable("u", SECOND);
    //    equation_systems("elastic").add_variable("v", SECOND);
    //    equation_systems("elastic").attach_assemble_function
(assemble_elastic);
    system.add_variable("u",FIRST);
    system.add_variable("v",FIRST);
    system.attach_assemble_function(assemble_elastic);

    // define equation_systems variable
    equation_systems.parameters.set<Real>("mu")=3D0.3;
    equation_systems.parameters.set<Real>("E")=3D20000;

    equation_systems.init();
    equation_systems.print_info();
    perf.stop_event("program init");

    perf.start_event("solve");
    equation_systems("elastic").solve();
    perf.stop_event("solve");

    perf.start_event("saving");
    GMVIO (mesh).write_equation_systems ("out.gmv", equation_systems);
    perf.stop_event("saving");
  }
  return libMesh::close();
}


void assemble_elastic(EquationSystems& es,
                      const std::string& system_name)
{
  assert (system_name =3D=3D "elastic");

  PerfLog perf("Matrix Assembly");
  const Mesh& mesh =3D es.get_mesh();
  //  const MeshData& mesh_data =3D es.get_mesh_data();
  const unsigned int dim =3D mesh.mesh_dimension();
  LinearImplicitSystem& system=3Des.get_system
<LinearImplicitSystem>("elastic");
  const Real mu =3D es.parameters.get<Real>("mu");
  const Real E =3D es.parameters.get<Real>("E");

  // Numeric ids corresponding to each variable in the system
  const unsigned int u_var=3Dsystem.variable_number("u");
  const unsigned int v_var=3Dsystem.variable_number("v");

  const DofMap& dof_map =3D system.get_dof_map();
  FEType fe_type =3D dof_map.variable_type(0);
  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
  // QGauss qrule (dim, FIFTH);
  QGauss qrule(dim,fe_type.default_quadrature_order());
  fe->attach_quadrature_rule (&qrule);
  // AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  // QGauss qface(dim-1, FIFTH);
  // fe_face->attach_quadrature_rule (&qface);
  const std::vector<Real>& JxW =3D fe->get_JxW();
  //  const std::vector<Point>& q_point =3D fe->get_xyz();
  //  const std::vector<std::vector<Real> >& phi =3D fe->get_phi();
  const std::vector<std::vector<RealGradient> >& dphi =3D fe->get_dphi();

  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  DenseSubMatrix<Number> Kuu(Ke),Kuv(Ke),
    Kvu(Ke),Kvv(Ke);
  DenseSubVector<Number> Fu(Fe),Fv(Fe);

  std::vector<unsigned int> dof_indices;
  std::vector<unsigned int> dof_indices_u;
  std::vector<unsigned int> dof_indices_v;

  MeshBase::const_element_iterator       el     =3D mesh.elements_begin();
  const MeshBase::const_element_iterator end_el =3D mesh.elements_end();
  for ( ; el !=3D end_el ; ++el){
    perf.start_event("elem init");
    const Elem* elem =3D *el;
    dof_map.dof_indices (elem, dof_indices);
    dof_map.dof_indices(elem,dof_indices_u,u_var);
    dof_map.dof_indices(elem,dof_indices_v,v_var);
    //
    const unsigned int n_dofs=3Ddof_indices.size();
    const unsigned int n_u_dofs=3Ddof_indices_u.size();
    const unsigned int n_v_dofs=3Ddof_indices_v.size();

    // Compute the element-specific data for the current element.
    fe->reinit (elem);
    Ke.resize (n_dofs,n_dofs);
    Fe.resize (n_dofs);
    // DenseSubMatrix.repositon()
    unsigned int ioff,joff;
    ioff=3D0;
    joff=3D0;
    Kuu.reposition(ioff,joff,n_u_dofs,n_u_dofs);
    joff +=3D n_u_dofs;
    Kuv.reposition(ioff,joff,n_u_dofs,n_v_dofs);
    ioff +=3D n_u_dofs;
    joff =3D 0;
    Kvu.reposition(ioff,joff,n_v_dofs,n_u_dofs);
    joff +=3D n_u_dofs;
    Kvv.reposition(ioff,joff,n_v_dofs,n_v_dofs);
    ioff=3D0;
    Fu.reposition(ioff,n_u_dofs);
    ioff +=3D n_u_dofs;
    Fv.reposition(ioff,n_v_dofs);

    perf.stop_event("elem init");

    perf.start_event("Ke");
    for (unsigned int qp=3D0; qp<qrule.n_points(); qp++){
      // Assemble Kuu
      for (unsigned int i=3D0; i<n_u_dofs;i++){
        for (unsigned int j=3D0; j<n_u_dofs;j++){
          Kuu(i,j) +=3D
JxW[qp]*(dphi[i][qp](0)*dphi[j][qp](0)/(1-mu*mu)+dphi[i][qp](1)*dphi[j][qp]=
(1)/2./(1.+mu))*E;
        }
      }
      // Assemble Kuv
      for(unsigned int i=3D0;i<n_u_dofs;i++){
        for(unsigned int j=3D0;j<n_v_dofs;j++){
          Kuv(i,j) +=3D
JxW[qp]*(dphi[i][qp](0)*dphi[j][qp](1)*mu/(1-mu*mu)+dphi[i][qp](1)*dphi[j][=
qp](0)/2./(1.+mu))*E;
        }
      }
      // Assemble Kvu
      for(unsigned int i=3D0;i<n_v_dofs;i++){
        for(unsigned int j=3D0;j<n_u_dofs;j++){
          Kvu(i,j) +=3D
JxW[qp]*(dphi[i][qp](1)*dphi[j][qp](0)*mu/(1-mu*mu)+dphi[i][qp](0)*dphi[j][=
qp](1)/2./(1.+mu))*E;
        }
      }
      // Assemble Kvv
      for(unsigned int i=3D0;i<n_v_dofs;i++){
        for(unsigned int j=3D0;j<n_v_dofs;j++){
          Kvv(i,j) +=3D
JxW[qp]*(dphi[i][qp](1)*dphi[j][qp](1)/(1-mu*mu)+dphi[i][qp](0)*dphi[j][qp]=
(0)/2./(1.+mu))*E;
        }
      }
    }
    perf.stop_event("Ke");
    // output Ke
    std::cout<<"output stiffness matrix"<<std::endl;
    for(unsigned int i=3D0;i<n_dofs;i++){
      for(unsigned int j=3D0;j<n_dofs;j++){
        std::cout<<"("<<i<<","<<j<<")"<<Ke(i,j)<<", ";
      }
      std::cout<<std::endl;
    }

    perf.start_event("matrix insertion");
    system.matrix->add_matrix (Ke, dof_indices);
    system.rhs->add_vector    (Fe, dof_indices);
    perf.stop_event("matrix insertion");
  }
  // assemble load and impose bc
  {
    assert(system_name =3D=3D "elastic");
    LinearImplicitSystem &system =3D es.get_system
<LinearImplicitSystem>("elastic");
    // Get writable references to the overall matrix and vector
    SparseMatrix<Number> &matrix =3D *system.matrix;
    NumericVector<Number> &rhs =3D *system.rhs;
    // Get a const reference to the mesh object
    const Mesh &mesh =3D es.get_mesh();
    // Apply load
    std::cout<<"Assemble Load"<<std::endl;
    unsigned int n_nodes =3D mesh.n_nodes();
    // fy=3D100 at node (1,0)
    for(unsigned int n_cnt=3D0;n_cnt<n_nodes;n_cnt++){
      const Node &curr_node=3Dmesh.node(n_cnt);
      const Real x=3Dcurr_node(0);
      const Real y=3Dcurr_node(1);
      /*
      if(fabs(y-0.0)<TOLERANCE && fabs(x-1.)<TOLERANCE){
        unsigned int ndof =3D curr_node.dof_number(0,1,0);
        rhs.set(ndof,100.0);  // fy=3D100
        break;
      }
      */
      //apply pressure=3D1 on x=3D1
      if(fabs(y-0.0)<TOLERANCE && fabs(x-1.)<TOLERANCE){
        unsigned int n =3D curr_node.dof_number(0,0,0);
        rhs.set(n,0.05);
      }
      if(fabs(y-1.)<TOLERANCE && fabs(x-1.)<TOLERANCE){
        unsigned int n =3D curr_node.dof_number(0,0,0);
        rhs.set(n,0.05);
      }
      if(0.01<y && y< 0.99 && fabs(x-1.)<TOLERANCE){
        unsigned int n =3D curr_node.dof_number(0,0,0);
        rhs.set(n,0.1);
      }
    }
    // fixed bc at left side
    for(unsigned int n_cnt=3D0;n_cnt<n_nodes;n_cnt++){
      const Node &curr_node=3Dmesh.node(n_cnt);
      const Real x=3Dcurr_node(0);
      const Real y=3Dcurr_node(1);
      const Real penalty =3D 1.e10;
      if(fabs(x-0.0)<TOLERANCE){
        unsigned int dofu =3D curr_node.dof_number(0,0,0);
        const Real u=3D0.0;
        rhs.set(dofu,u*penalty);

std::cout<<"dofu=3D"<<dofu<<std::endl<<"matrix"<<matrix(dofu,dofu)<<std::en=
dl;

        //        matrix.set(dofu,dofu,penalty);
        matrix.add(dofu,dofu,penalty);
        /*
        unsigned int dofv =3D curr_node.dof_number(0,1,0);
        const Real v=3D0.0;
        rhs.set(dofv,v*penalty);
        matrix.add(dofv,dofv,penalty);
        */
        if(fabs(y-0.0)<TOLERANCE){
          unsigned int n =3D curr_node.dof_number(0,1,0);
          const Real v=3D0.0;
          rhs.set(n,v*penalty);
          matrix.add(n,n,penalty);
        }
      }
    }
  }
  {
    apply_periodic_bc(es,system_name);
  }
}

void apply_periodic_bc(EquationSystems &es, const std::string &system_name)=
{
  assert(system_name=3D=3D"elastic");
  PerfLog perf_log("Periodic bc. Assembly",false);
  const Mesh &mesh=3Des.get_mesh();
  LinearImplicitSystem &system=3Des.get_system
<LinearImplicitSystem>("elastic");
  DofMap &dof_map=3Dsystem.get_dof_map();

  // impose ux(0.5,0.0)=3Dux(0.8,0.0)
  for(unsigned int n=3D0;n<mesh.n_nodes();n++){
    const Node &node1=3Dmesh.node(n);
    const double x1=3Dnode1(0);
    if((fabs(node1(0)-0.5)<TOLERANCE) && fabs(node1(1))<TOLERANCE){
      for(unsigned int n1=3D0;n1<mesh.n_nodes();n1++){
        const Node &node2=3Dmesh.node(n1);
        const double x2=3Dnode2(0);
        std::cout<<"x1=3D"<<x1<<",x2=3D"<<x2<<std::endl;
        if(fabs(node2(0)-0.8)<TOLERANCE && fabs(node2(1))<TOLERANCE){
          unsigned int d1=3Dnode1.dof_number(0,0,0);
          unsigned int d2=3Dnode2.dof_number(0,0,0);
          std::cout<<"n_left=3D"<<n<<",n_right=3D"<<n1<<std::endl;
          DofConstraintRow constraint;
          constraint[d1]=3D1;
          dof_map.add_constraint_row(d2,constraint);
          break;
        }
      }
    }
  }
}