Hello Daniel,
Thanks for the reply.
I am actually doing what you suggested. Please look at the updated code
that I had posted earlier and works with deal.ii 8.5.
These lines are already there in the code. Surprisingly, such an error is
still showing up. I am not sure what is still causing this weird behavior
(bug?).
IndexSet locally_relevant_dofs2 = locally_relevant_dofs;
locally_relevant_dofs2.add_index(the_dof);
constraint_temp.reinit(locally_relevant_dofs2);
I am attaching the code again for your perusal.
On Thursday, April 20, 2017 at 5:37:46 AM UTC-4, Daniel Arndt wrote:
>
> Rajat,
>
> What is happening here is that you try to create constraints that include
> DoFs that are not element of the IndexSet you provided to the
> DynamicSparsityPattern object.
> Therefore, all the entries for your "fixed" DoF are ignored if it is not
> already part of the locally relevant DoFs.
> Have a look at the discussion in "The LaplaceProblem class;
> LaplaceProblem::setup_system" in step-40 [1].
>
> You can solve your issue by changing
> IndexSet locally_relevant_dofs2 = locally_relevant_dofs;
> constraint_temp.reinit(locally_relevant_dofs2);
> to
> IndexSet locally_relevant_dofs2 = locally_relevant_dofs;
> locally_relevant_dofs2.add_index(the_dof);
> constraint_temp.reinit(locally_relevant_dofs2);
>
> Best,
> Daniel
>
> [1] http://www.dealii.org/8.4.1/doxygen/deal.II/step_40.html
>
--
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#ifndef COMMON_H
#define COMMON_H
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/generic_linear_algebra.h>
#define USE_PETSC_LA
namespace LA
{
#ifdef USE_PETSC_LA
using namespace dealii::LinearAlgebraPETSc;
#else
using namespace dealii::LinearAlgebraTrilinos;
#endif
}
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/distributed/solution_transfer.h>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <dirent.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
#include <unistd.h>
#endif
#include "common.h"
using namespace dealii;
template<int dim>
class BugFinder
{
private:
parallel::distributed::Triangulation<dim> triangulation;
FESystem<dim> fe; // 3 dofs per node
DoFHandler<dim> dof_handler;
IndexSet locally_relevant_dofs;
IndexSet locally_owned_dofs;
LA::MPI::Vector system_rhs;
LA::MPI::SparseMatrix system_matrix;
ConstraintMatrix constraints;
std::vector<unsigned int> b_c_codes;
public:
BugFinder();
~BugFinder() {
dof_handler.clear();
};
void apply_constraints();
void assemble_system();
void output_data();
};
template<int dim>
void BugFinder<dim>::assemble_system()
{
unsigned int dofs_per_cell = fe.dofs_per_cell;
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<unsigned int> local_dof_indices(dofs_per_cell);
system_matrix = 0;
system_rhs = 0;
for (; cell != endc; ++cell)
if (cell->is_locally_owned())
{
for (int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs[i] = std::rand();
for (int j = 0; j < dofs_per_cell; ++j)
{
cell_matrix[i][j] = std::rand();
}
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
system_matrix.compress (VectorOperation::add);
}
template<int dim>
void BugFinder<dim>::output_data()
{
std::cout << "output_data ... " ;
std::string name = "grid";
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
Vector<float> subdomain (triangulation.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
{
subdomain(i) = triangulation.locally_owned_subdomain();
//std::cout<<"subdomain("<<i<<") ="<< subdomain(i)<<std::endl;
}
data_out.add_data_vector (subdomain, "subdomain");
data_out.build_patches ();
const std::string filename = (name +
"." +
Utilities::int_to_string
(triangulation.locally_owned_subdomain()));
std::ofstream output ((filename + ".vtu").c_str());
data_out.write_vtu (output);
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{
std::vector<std::string> filenames;
for (unsigned int i = 0;
i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
++i)
filenames.push_back (name +
"." +
Utilities::int_to_string (i) +
".vtu");
std::ofstream master_output ((filename + ".pvtu").c_str());
data_out.write_pvtu_record (master_output, filenames);
}
std::cout << " ends" << std::endl;
}
template<int dim>
void BugFinder<dim>::apply_constraints() // fix all the z components on the +z face to a single component given by the dof
{
ConstraintMatrix constraint_temp;
constraints.clear(); constraint_temp.clear();
constraints.reinit(locally_relevant_dofs);
constraint_temp.reinit(locally_relevant_dofs);
unsigned int the_dof = 0;
unsigned int plus_z_face = b_c_codes[2];
typename DoFHandler<dim>::active_cell_iterator
v_cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; v_cell != endc; ++v_cell)
if (v_cell->is_locally_owned())
if (v_cell->at_boundary())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
if (v_cell->face(face)->at_boundary() && v_cell->face(face)->boundary_id() == plus_z_face)
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
{
if (locally_owned_dofs.is_element(v_cell->face(face)->vertex_dof_index(v, 2)))
{
the_dof = std::max(v_cell->face(face)->vertex_dof_index(v, 2), the_dof);
}
}
the_dof = Utilities::MPI::max(the_dof, MPI_COMM_WORLD);
ComponentMask cm(dim, false);
cm.set(2, true); // z component
IndexSet selected_dofs;
std::set< types::boundary_id > boundary_ids;
boundary_ids.insert(plus_z_face);
DoFTools::extract_boundary_dofs (dof_handler, cm, selected_dofs, boundary_ids);
std::vector<types::global_dof_index> index_vector;
// selected_dofs.print(std::cout);
selected_dofs.fill_index_vector(index_vector);
IndexSet locally_relevant_dofs2 = locally_relevant_dofs;
locally_relevant_dofs2.add_index(the_dof);
constraint_temp.reinit(locally_relevant_dofs2);
constraints.reinit(locally_relevant_dofs2);
constraint_temp.get_local_lines().print(std::cout ); MPI_Barrier(MPI_COMM_WORLD);
for (unsigned int i = 0; i < index_vector.size(); ++i)
{
if (index_vector[i] == the_dof)
continue;
if (!locally_relevant_dofs.is_element(index_vector[i]))
continue;
std::cout << " dof constrained " << index_vector[i]
//<< " lrd ? " << locally_relevant_dofs.is_element(index_vector[i])
<< std::endl;
constraint_temp.add_line(index_vector[i]);
constraint_temp.add_entry(index_vector[i], the_dof, 1.0);
}
constraint_temp.close();
constraints.merge(constraint_temp);
constraints.close ();
MPI_Barrier(MPI_COMM_WORLD);;
constraints.print(std::cout);
MPI_Barrier(MPI_COMM_WORLD);;
std::cout << " constraints printed " <<std::endl;
DynamicSparsityPattern dsp (locally_relevant_dofs);
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
SparsityTools::distribute_sparsity_pattern (dsp,
dof_handler.n_locally_owned_dofs_per_processor(),
MPI_COMM_WORLD,
locally_relevant_dofs);
system_matrix.reinit (locally_owned_dofs,
locally_owned_dofs,
dsp,
MPI_COMM_WORLD);
system_rhs.reinit(locally_owned_dofs, MPI_COMM_WORLD);
}
template<int dim>
BugFinder<dim>::BugFinder():
triangulation (MPI_COMM_WORLD,
typename Triangulation<dim>::MeshSmoothing
(Triangulation<dim>::smoothing_on_refinement
| Triangulation<dim>::smoothing_on_coarsening),
parallel::distributed::Triangulation<dim>::no_automatic_repartitioning),
fe(FE_Q<dim>(1), dim),
dof_handler(triangulation)
{
double Point1[] = {.1, .1, 0.3};
double Point2[] = { -.1, -.1, 0};
double limits[3][2];
Point<dim> p1(Point1[0], Point1[1], Point1[2]);
Point<dim> p2(Point2[0], Point2[1], Point2[2]);
for (int d = 0; d < dim; ++d)
{
limits[d][0] = p1[d] > p2[d] ? p2[d] : p1[d];
limits[d][1] = p1[d] > p2[d] ? p1[d] : p2[d];
}
std::vector<unsigned int> repetitions;
repetitions.push_back(4);
repetitions.push_back(4);
repetitions.push_back(4);
GridGenerator::subdivided_hyper_rectangle(triangulation,
repetitions,
p1, p2);
dof_handler.distribute_dofs(fe);
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
b_c_codes.push_back(101);
b_c_codes.push_back(2);
b_c_codes.push_back(3);
/*****************************************************/
{ typename Triangulation<dim>::cell_iterator
cell = triangulation.begin(),
endc = triangulation.end();
for (; cell != endc; ++cell)
{
for (unsigned int face = 0;
face < GeometryInfo<dim>::faces_per_cell; ++face)
{
if(cell->face(face)->at_boundary())
{
double negative_limit_this_dimension = limits[2][0];
double positive_limit_this_dimension = limits[2][1];
// std::cout << "negative_limit_this_dimension " << negative_limit_this_dimension << " d " << d << std::endl;
// std::cout << "positive_limit_this_dimension " << positive_limit_this_dimension << " d " << d << std::endl;
if (std::fabs(cell->face(face)->center()(2)
- negative_limit_this_dimension) < 1e-12)
cell->face(face)->set_boundary_id(b_c_codes[1]);
else if (std::fabs(cell->face(face)->center()(2)
- positive_limit_this_dimension) < 1e-12)
cell->face(face)->set_boundary_id(b_c_codes[2]);
else
cell->face(face)->set_boundary_id(b_c_codes[0]);
}
}
}
}
/*****************************************************/
apply_constraints(); std::cout << "constraints done " << std::endl;
assemble_system();
}
int main(int argc, char *argv[])
{
dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
BugFinder<3> bugfinder;
std::cout << "created " << std::endl;
return 0;
}
##
# CMake script for the fdm tutorial program:
##
# Set the name of the project and target:
SET(TARGET "fdm")
add_compile_options(-Wno-comment -Wno-unused-parameter -Wno-unused)
# Declare all source files the target consists of:
include_directories("${CMAKE_SOURCE_DIR}/src")
#FILE(GLOB SRC_FILES ./src/main.cc ./src/velocity_equation.cc)
SET(TARGET_SRC ./src/main.cc)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.0 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
#
# Are all dependencies fullfilled?
#
IF(NOT DEAL_II_WITH_PETSC OR NOT DEAL_II_WITH_P4EST)
MESSAGE(FATAL_ERROR "
Error! The deal.II library found at ${DEAL_II_PATH} was not configured with
DEAL_II_WITH_PETSC = ON
DEAL_II_WITH_P4EST = ON
One or all of these are OFF in your installation but are required for this
tutorial step."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET} CXX)
DEAL_II_INVOKE_AUTOPILOT()
