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=u2. 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=0, p=1 on
x=1, 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  USA */

#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=1; 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 =
      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")=0.3;
    equation_systems.parameters.set<Real>("E")=20000;

    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 == "elastic");

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

  // Numeric ids corresponding to each variable in the system
  const unsigned int u_var=system.variable_number("u");
  const unsigned int v_var=system.variable_number("v");
 
  const DofMap& dof_map = system.get_dof_map();
  FEType fe_type = 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 = fe->get_JxW();
  //  const std::vector<Point>& q_point = fe->get_xyz();
  //  const std::vector<std::vector<Real> >& phi = fe->get_phi();
  const std::vector<std::vector<RealGradient> >& dphi = 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     = mesh.elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.elements_end();
  for ( ; el != end_el ; ++el){
    perf.start_event("elem init");
    const Elem* elem = *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=dof_indices.size();
    const unsigned int n_u_dofs=dof_indices_u.size();
    const unsigned int n_v_dofs=dof_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=0;
    joff=0;
    Kuu.reposition(ioff,joff,n_u_dofs,n_u_dofs);
    joff += n_u_dofs;
    Kuv.reposition(ioff,joff,n_u_dofs,n_v_dofs);
    ioff += n_u_dofs;
    joff = 0;
    Kvu.reposition (ioff,joff,n_v_dofs,n_u_dofs);
    joff += n_u_dofs;
    Kvv.reposition(ioff,joff,n_v_dofs,n_v_dofs);
    ioff=0;
    Fu.reposition(ioff,n_u_dofs);
    ioff += n_u_dofs;
    Fv.reposition(ioff,n_v_dofs);

    perf.stop_event("elem init");

    perf.start_event("Ke");
    for (unsigned int qp=0; qp<qrule.n_points(); qp++){
      // Assemble Kuu
      for (unsigned int i=0; i<n_u_dofs;i++){
        for (unsigned int j=0; j<n_u_dofs;j++){
          Kuu(i,j) += 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=0;i<n_u_dofs;i++){
        for(unsigned int j=0;j<n_v_dofs;j++){
          Kuv(i,j) += 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=0;i<n_v_dofs;i++){
        for(unsigned int j=0;j<n_u_dofs;j++){
          Kvu(i,j) += 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=0;i<n_v_dofs;i++){
        for(unsigned int j=0;j<n_v_dofs;j++){
          Kvv(i,j) += 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=0;i<n_dofs;i++){
      for(unsigned int j=0;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 == "elastic");
    LinearImplicitSystem &system = es.get_system<LinearImplicitSystem>("elastic");
    // Get writable references to the overall matrix and vector
    SparseMatrix<Number> &matrix = *system.matrix;
    NumericVector<Number> &rhs = *system.rhs;
    // Get a const reference to the mesh object
    const Mesh &mesh = es.get_mesh();
    // Apply load
    std::cout<<"Assemble Load"<<std::endl;
    unsigned int n_nodes = mesh.n_nodes();
    // fy=100 at node (1,0)
    for(unsigned int n_cnt=0;n_cnt<n_nodes;n_cnt++){
      const Node &curr_node=mesh.node(n_cnt);
      const Real x=curr_node(0);
      const Real y=curr_node(1);
      /*
      if(fabs(y-0.0)<TOLERANCE && fabs(x-1.)<TOLERANCE){
        unsigned int ndof = curr_node.dof_number(0,1,0);
        rhs.set(ndof,100.0);  // fy=100
        break;
      }
      */
      //apply pressure=1 on x=1
      if(fabs(y-0.0)<TOLERANCE && fabs(x-1.)<TOLERANCE){
        unsigned int n = curr_node.dof_number(0,0,0);
        rhs.set(n,0.05);
      }
      if(fabs(y-1.)<TOLERANCE && fabs(x-1.)<TOLERANCE){
        unsigned int n = 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 = curr_node.dof_number(0,0,0);
        rhs.set(n,0.1);
      }
    }
    // fixed bc at left side
    for(unsigned int n_cnt=0;n_cnt<n_nodes;n_cnt++){
      const Node &curr_node=mesh.node(n_cnt);
      const Real x=curr_node(0);
      const Real y=curr_node(1);
      const Real penalty = 1.e10;
      if(fabs(x-0.0)<TOLERANCE){
        unsigned int dofu = curr_node.dof_number(0,0,0);
        const Real u=0.0;
        rhs.set(dofu,u*penalty);
        std::cout<<"dofu="<<dofu<<std::endl<<"matrix"<<matrix(dofu,dofu)<<std::endl;

        //        matrix.set(dofu,dofu,penalty);
        matrix.add(dofu,dofu,penalty);
        /*
        unsigned int dofv = curr_node.dof_number(0,1,0);
        const Real v=0.0;
        rhs.set(dofv,v*penalty);
        matrix.add(dofv,dofv,penalty);
        */
        if(fabs(y-0.0)<TOLERANCE){
          unsigned int n = curr_node.dof_number(0,1,0);
          const Real v=0.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=="elastic");
  PerfLog perf_log("Periodic bc. Assembly",false);
  const Mesh &mesh= es.get_mesh();
  LinearImplicitSystem &system=es.get_system<LinearImplicitSystem>("elastic");
  DofMap &dof_map=system.get_dof_map();

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