From 8e68a19ec05d6b6241eef46d95853e6a7cce7d74 Mon Sep 17 00:00:00 2001 From: "denis.davydov" Date: Tue, 20 May 2014 19:58:13 +0000 Subject: [PATCH] added a test of matrix_scalar_product git-svn-id: https://svn.dealii.org/branches/branch_petscscalar_complex@32949 0785d39b-7218-0410-832d-ea1e28bc413d --- .../parallel_sparse_matrix_01.cc | 210 ++++++++++++++++++ .../parallel_sparse_matrix_01.output | 5 + 2 files changed, 215 insertions(+) create mode 100644 tests/petsc_complex/parallel_sparse_matrix_01.cc create mode 100644 tests/petsc_complex/parallel_sparse_matrix_01.output diff --git a/tests/petsc_complex/parallel_sparse_matrix_01.cc b/tests/petsc_complex/parallel_sparse_matrix_01.cc new file mode 100644 index 0000000000..c184f0871e --- /dev/null +++ b/tests/petsc_complex/parallel_sparse_matrix_01.cc @@ -0,0 +1,210 @@ +// --------------------------------------------------------------------- +// $Id: testsuite.html 32446 2014-02-10 14:31:22Z bangerth $ +// +// Copyright (C) 2013 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// test that matrix_scalar_product of a symmetric matrix +// applied to the same vector result in a real number + + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +using namespace dealii; + +static const unsigned int dim = 2; + +void test () +{ + + MPI_Comm mpi_communicator (MPI_COMM_WORLD); + const unsigned int n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)); + const unsigned int this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)); + + ConditionalOStream pcout (std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) + == 0)); + + parallel::distributed::Triangulation tria(mpi_communicator, + typename Triangulation::MeshSmoothing + (Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)); + + if (false) + { + Triangulation triangulation1; + Triangulation triangulation2; + Triangulation triangulation12; + GridGenerator::hyper_cube (triangulation1, -1,0); //create a square [-1,0]^d domain + GridGenerator::hyper_cube (triangulation2, -1,0); //create a square [-1,0]^d domain + Point shift_vector; + shift_vector[0] = 1.0; + GridTools::shift(shift_vector,triangulation2); + GridGenerator::merge_triangulations (triangulation1, triangulation2, tria); + } + else + { + GridGenerator::hyper_cube (tria, -1,0); + //tria.refine_global (1); + } + + const unsigned int poly_degree = 1; + FE_Q fe(poly_degree); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs (); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs); + + PETScWrappers::MPI::Vector vector, conjugate; + PETScWrappers::MPI::SparseMatrix mass_matrix; + + + vector.reinit(locally_owned_dofs, mpi_communicator); + const types::global_dof_index n_local_dofs = locally_owned_dofs.n_elements(); + mass_matrix.reinit (mpi_communicator, + dof_handler.n_dofs(), + dof_handler.n_dofs(), + n_local_dofs, + n_local_dofs, + dof_handler.max_couplings_between_dofs()); + + ConstraintMatrix constraints; + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + //assemble mass matrix: + mass_matrix = 0.0; + { + //a separate quadrature formula + //enough for mass and kinetic matrices assembly + QGauss quadrature_formula(poly_degree+1); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + cell_mass_matrix = PetscScalar(0); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints + .distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } + mass_matrix.compress (VectorOperation::add); + } + + vector = std::complex(1.0,-1.0); + + //deallog<<"original vector: "< norm = + mass_matrix.matrix_scalar_product(vector, vector); + deallog<