Hello All !

I was trying the step-3 
<https://www.dealii.org/current/doxygen/deal.II/step_3.html> with periodic 
boundary conditions and dim = 2 and spacedim = 3. The main reason behind 
this is I want to test periodic boundary conditions for dim != spacedim. I 
have tested the code for dim = spacedim = 2 and the periodic boundary 
conditions and it works (never mind the results for now !). I have also 
tested the code for dim = 2 and spacedim = 3 without periodic boundary 
conditions and it works too. But only when I add periodic conditions for 
that case, the issue arises. 

I got the following error,

error: undefined reference to 'void 
dealii::DoFTools::make_periodicity_constraints<2, 3, 
double>(dealii::DoFHandler<2, 3> const&, unsigned int, unsigned int, unsigned 
int, dealii::AffineConstraints<double>&, dealii::ComponentMask const&, double)'
collect2: error: ld returned 1 exit status
make[2]: *** [CMakeFiles/step1.dir/build.make:139: step1] Error 1
make[1]: *** [CMakeFiles/Makefile2:76: CMakeFiles/step1.dir/all] Error 2
make: *** [Makefile:84: all] Error 2

The version that I am using is 9.3.0-pre. Although I am not really great at 
C++ but according to the template 
<https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/include/deal.II/dofs/dof_tools.h#L1144>
 
and source 
<https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/source/dofs/dof_tools_constraints.cc#L2552>,
 
it should have worked, but I could not find anything here 
<https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/source/dofs/dof_tools_constraints.inst.in#L69>.
 
I tried running it without template arguments or with only '2' as a 
template argument, but couldn't make it work. I have attached the code 
which I have tried, and also the code dim = 2 and spacedim = 3 without 
periodic boundary conditions. Any help will be greatly appreciated.

Thank You

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to [email protected].
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e11a9f41-4750-459d-bedc-ef0b3ae83dddn%40googlegroups.com.
/* ---------------------------------------------------------------------
 *
 * Copyright (C) 1999 - 2020 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, 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.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * ---------------------------------------------------------------------
 *
 * Authors: Wolfgang Bangerth, 1999,
 *          Guido Kanschat, 2011
 */
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
using namespace dealii;
class Step3
{
public:
  Step3();
  void run();
private:
  void make_grid();
  void setup_system();
  void assemble_system();
  void solve();
  void output_results() const;
  Triangulation<2,3> triangulation;
  FE_Q<2,3>          fe;
  DoFHandler<2,3>    dof_handler;
  MappingQ<2,3>      mapping;
  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;
  Vector<double> solution;
  Vector<double> system_rhs;
};
Step3::Step3()
  : fe(1)
  , dof_handler(triangulation)
  , mapping(2)
{}
void Step3::make_grid()
{
  GridGenerator::subdivided_hyper_cube(triangulation, 10, -10.0, 10.0, 1);
  // triangulation.refine_global(5);
  std::cout << "Number of active cells: " << triangulation.n_active_cells()
            << std::endl;
}
void Step3::setup_system()
{
  dof_handler.distribute_dofs(fe);
  std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
            << std::endl;
  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler, dsp);
  sparsity_pattern.copy_from(dsp);
  system_matrix.reinit(sparsity_pattern);
  solution.reinit(dof_handler.n_dofs());
  system_rhs.reinit(dof_handler.n_dofs());
}
void Step3::assemble_system()
{
  QGauss<2> quadrature_formula(fe.degree + 1);
  FEValues<2,3> fe_values(mapping, fe,
                        quadrature_formula,
                        update_values | update_gradients | update_JxW_values);
  const unsigned int dofs_per_cell = fe.dofs_per_cell;
  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
  Vector<double>     cell_rhs(dofs_per_cell);
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
  for (const auto &cell : dof_handler.active_cell_iterators())
    {
      fe_values.reinit(cell);
      cell_matrix = 0;
      cell_rhs    = 0;
      for (const unsigned int q_index : fe_values.quadrature_point_indices())
        {
          for (const unsigned int i : fe_values.dof_indices())
            for (const unsigned int j : fe_values.dof_indices())
              cell_matrix(i, j) +=
                (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
                 fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
                 fe_values.JxW(q_index));           // dx
          for (const unsigned int i : fe_values.dof_indices())
            cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
                            1 *                                 // f(x_q)
                            fe_values.JxW(q_index));            // dx
        }
      cell->get_dof_indices(local_dof_indices);
      for (const unsigned int i : fe_values.dof_indices())
        for (const unsigned int j : fe_values.dof_indices())
          system_matrix.add(local_dof_indices[i],
                            local_dof_indices[j],
                            cell_matrix(i, j));
      for (const unsigned int i : fe_values.dof_indices())
        system_rhs(local_dof_indices[i]) += cell_rhs(i);
    }
  std::map<types::global_dof_index, double> boundary_values;
  VectorTools::interpolate_boundary_values(mapping, dof_handler,
                                           0,
                                           Functions::ZeroFunction<3>(),
                                           boundary_values);
  MatrixTools::apply_boundary_values(boundary_values,
                                     system_matrix,
                                     solution,
                                     system_rhs);
}
void Step3::solve()
{
  SolverControl solver_control(1000, 1e-12);
  SolverCG<Vector<double>> solver(solver_control);
  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
void Step3::output_results() const
{
  DataOut<2, DoFHandler<2,3>> data_out;
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(solution, "solution");
  data_out.build_patches();
  std::ofstream output("solution.vtk");
  data_out.write_vtk(output);
}
void Step3::run()
{
  make_grid();
  setup_system();

  {
    FEValuesExtractors::Scalar u(0);
    std::vector<GridTools::PeriodicFacePair<typename DoFHandler<2,3>::cell_iterator>> periodicity_vector;
    const unsigned int direction = 0;
    GridTools::collect_periodic_faces(dof_handler, 0, 1, direction, periodicity_vector);
    std::cout << "\nSize of periodicity_vector : " << periodicity_vector.size() << std::endl;
    AffineConstraints<double> constraints;
    constraints.clear();
    DoFTools::make_periodicity_constraints<2,3>(dof_handler, 0, 1, direction, constraints, fe.component_mask(u));
    constraints.close();
  }

  assemble_system();
  solve();
  output_results();
}
int main()
{
  deallog.depth_console(2);
  Step3 laplace_problem;
  laplace_problem.run();
  return 0;
}
/* ---------------------------------------------------------------------
 *
 * Copyright (C) 1999 - 2020 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, 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.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * ---------------------------------------------------------------------
 *
 * Authors: Wolfgang Bangerth, 1999,
 *          Guido Kanschat, 2011
 */
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
using namespace dealii;
class Step3
{
public:
  Step3();
  void run();
private:
  void make_grid();
  void setup_system();
  void assemble_system();
  void solve();
  void output_results() const;
  Triangulation<2,3> triangulation;
  FE_Q<2,3>          fe;
  DoFHandler<2,3>    dof_handler;
  MappingQ<2,3>      mapping;
  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;
  Vector<double> solution;
  Vector<double> system_rhs;
};
Step3::Step3()
  : fe(1)
  , dof_handler(triangulation)
  , mapping(2)
{}
void Step3::make_grid()
{
  GridGenerator::subdivided_hyper_cube(triangulation, 10, -10.0, 10.0, 1);
  // triangulation.refine_global(5);
  std::cout << "Number of active cells: " << triangulation.n_active_cells()
            << std::endl;
}
void Step3::setup_system()
{
  dof_handler.distribute_dofs(fe);
  std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
            << std::endl;
  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler, dsp);
  sparsity_pattern.copy_from(dsp);
  system_matrix.reinit(sparsity_pattern);
  solution.reinit(dof_handler.n_dofs());
  system_rhs.reinit(dof_handler.n_dofs());
}
void Step3::assemble_system()
{
  QGauss<2> quadrature_formula(fe.degree + 1);
  FEValues<2,3> fe_values(mapping, fe,
                        quadrature_formula,
                        update_values | update_gradients | update_JxW_values);
  const unsigned int dofs_per_cell = fe.dofs_per_cell;
  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
  Vector<double>     cell_rhs(dofs_per_cell);
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
  for (const auto &cell : dof_handler.active_cell_iterators())
    {
      fe_values.reinit(cell);
      cell_matrix = 0;
      cell_rhs    = 0;
      for (const unsigned int q_index : fe_values.quadrature_point_indices())
        {
          for (const unsigned int i : fe_values.dof_indices())
            for (const unsigned int j : fe_values.dof_indices())
              cell_matrix(i, j) +=
                (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
                 fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
                 fe_values.JxW(q_index));           // dx
          for (const unsigned int i : fe_values.dof_indices())
            cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
                            1 *                                 // f(x_q)
                            fe_values.JxW(q_index));            // dx
        }
      cell->get_dof_indices(local_dof_indices);
      for (const unsigned int i : fe_values.dof_indices())
        for (const unsigned int j : fe_values.dof_indices())
          system_matrix.add(local_dof_indices[i],
                            local_dof_indices[j],
                            cell_matrix(i, j));
      for (const unsigned int i : fe_values.dof_indices())
        system_rhs(local_dof_indices[i]) += cell_rhs(i);
    }
  std::map<types::global_dof_index, double> boundary_values;
  VectorTools::interpolate_boundary_values(mapping, dof_handler,
                                           0,
                                           Functions::ZeroFunction<3>(),
                                           boundary_values);
  MatrixTools::apply_boundary_values(boundary_values,
                                     system_matrix,
                                     solution,
                                     system_rhs);
}
void Step3::solve()
{
  SolverControl solver_control(1000, 1e-12);
  SolverCG<Vector<double>> solver(solver_control);
  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
void Step3::output_results() const
{
  DataOut<2, DoFHandler<2,3>> data_out;
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(solution, "solution");
  data_out.build_patches();
  std::ofstream output("solution.vtk");
  data_out.write_vtk(output);
}
void Step3::run()
{
  make_grid();
  setup_system();

  {
    FEValuesExtractors::Scalar u(0);
    std::vector<GridTools::PeriodicFacePair<typename DoFHandler<2,3>::cell_iterator>> periodicity_vector;
    const unsigned int direction = 0;
    GridTools::collect_periodic_faces(dof_handler, 0, 1, direction, periodicity_vector);
    std::cout << "\nSize of periodicity_vector : " << periodicity_vector.size() << std::endl;
    // AffineConstraints<double> constraints;
    // constraints.clear();
    // DoFTools::make_periodicity_constraints<2,3>(dof_handler, 0, 1, direction, constraints, fe.component_mask(u));
    // constraints.close();
  }

  assemble_system();
  solve();
  output_results();
}
int main()
{
  deallog.depth_console(2);
  Step3 laplace_problem;
  laplace_problem.run();
  return 0;
}

Reply via email to