Hi all,
I wrote a program where I use a HyperShell.
Next, I produce the right-hand side of my system by computing a surface
integral. In this Integral, I integrate a constant function times a basis
function.
One expects the resulting vector to be constant but I get different values
in it.
I have visualized the right-hand side and have attached the figure here.
I have also made and attached a minimal program that produces this figure.
I don't think this behavior is correct. Does anyone have any ideas? Have I
missed anything in the implementation or otherwise?
Best,
Ashkan
--
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.
# top level CMakeLists.txt
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.0 REQUIRED
HINTS
${DEAL_II_DIR} $ENV{DEAL_II_DIR} )
DEAL_II_INITIALIZE_CACHED_VARIABLES()
# Set project name and the languages used in it
PROJECT(example3DSphere CXX)
########################## V A R I A B L E S ##########################
# Set the executable name
SET(EXEC_NAME ex)
# set which files to clean in the clear target
SET(CLEAN_UP_FILES
*gmv *gnuplot *gpl *eps *pov *vtk *vtu *ucd *.d2 *dat *.log *.m *.1 *.mod
*.log
)
ADD_EXECUTABLE(${EXEC_NAME} main.cc)
########################### T A R G E T S #############################
ADD_CUSTOM_TARGET(clear
COMMAND rm -f ${CLEAN_UP_FILES}
WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
COMMENT "Cleanup output created by running the application"
)
ADD_CUSTOM_TARGET(debug
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR}
COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
COMMENT "Switch to Debug mode"
)
ADD_CUSTOM_TARGET(release
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
COMMENT "Switch to Release mode"
)
DEAL_II_SETUP_TARGET(${EXEC_NAME})
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/std_cxx11/bind.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/utilities.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/dofs/block_info.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/data_out_faces.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/index_set.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include "deal.II/base/utilities.h"
#include "deal.II/base/logstream.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <typeinfo>
#include <vector>
using namespace dealii;
using std::endl;
using std::stringstream;
using std::string;
using std::ostringstream;
using std::ofstream;
namespace LA {
using namespace dealii::LinearAlgebraTrilinos;
}
#define DIM 3
namespace params {
double earth_inner = 3.481e6;
double earth_outer = 6.371e6;
int n_refs = 4;
int degree = 1;
}
class spGeom{
public:
spGeom();
~spGeom();
void run();
private:
MPI_Comm mpi_communicator;
ConditionalOStream pcout;
TimerOutput timer;
parallel::distributed::Triangulation<DIM> triangulation;
DoFHandler<DIM> dof_handler;
FESystem<DIM> fe;
ConstraintMatrix constraints;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
std::vector<IndexSet> owned_partitioning;
std::vector<IndexSet> relevant_partitioning;
LA::MPI::BlockSparseMatrix system_matrix;
LA::MPI::BlockVector system_rhs;
std::vector<types::global_dof_index> dofs_per_block;
// Methods
void makeGrid();
void assemble();
void setupDofs();
void output_results();
};
spGeom::spGeom() :
mpi_communicator(MPI_COMM_WORLD),
pcout(std::cout,
(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
timer(mpi_communicator,
pcout,
TimerOutput::summary,
TimerOutput::wall_times),
triangulation (mpi_communicator,
typename Triangulation<DIM>::MeshSmoothing
(Triangulation<DIM>::smoothing_on_refinement |
Triangulation<DIM>::smoothing_on_coarsening)),
dof_handler (triangulation),
fe (FE_Q<DIM>(params::degree+1), DIM, // Displacement
FE_Q<DIM>(params::degree), 1), // pressure
dofs_per_block(std::vector<types::global_dof_index>(2))
{}
spGeom::~spGeom()
{
dof_handler.clear();
}
void spGeom::makeGrid(){
pcout << "Make grid" << endl;
TimerOutput::Scope t(timer, "Grid generator");
// Number of initial subdivisions for each axis
static SphericalManifold<DIM> surface_description;
GridGenerator::hyper_shell(triangulation,
Point<DIM>(), // center is at zero
params::earth_inner,
params::earth_outer,
12,
true);
triangulation.set_all_manifold_ids(0);
triangulation.set_manifold(0, surface_description);
// Refine the mesh with the number of refinement in the parameters.
triangulation.refine_global(params::n_refs);
}
void spGeom::setupDofs(){
pcout << "Setup" << endl;
// Set scope of the timer
TimerOutput::Scope t(timer, "Setup DOFs");
dof_handler.distribute_dofs (fe);
// Renumber component wise
std::vector<unsigned int> gcomponents (DIM+1, 0);
gcomponents[DIM] = 1;
// DOF renumbering
DoFRenumbering::component_wise (dof_handler, gcomponents);
// Get number of dofs in components
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, gcomponents);
// Output number of dofs
const unsigned int n_u = dofs_per_block[0],
n_p = dofs_per_block[1];
locally_owned_dofs = dof_handler.locally_owned_dofs();
owned_partitioning.clear();
relevant_partitioning.clear();
owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u+n_p));
DoFTools::extract_locally_relevant_dofs (dof_handler,
locally_relevant_dofs);
relevant_partitioning.push_back(locally_relevant_dofs.get_view(0,n_u));
relevant_partitioning.push_back(locally_relevant_dofs.get_view(n_u,n_u+n_p));
//Interpolate boudaries using constraint matrix
{
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
ComponentMask ns_mask = fe.component_mask(FEValuesExtractors::Vector(0)); //NO_SLIP
// set inner boundary to non moving
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction<DIM>(DIM+1),
constraints,
ns_mask);
constraints.close();
}
system_matrix.clear();
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints, true);
SparsityTools::distribute_sparsity_pattern (dsp,
dof_handler.locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit(owned_partitioning,
dsp,
mpi_communicator);
system_rhs.reinit (owned_partitioning,
mpi_communicator);
}
void spGeom::assemble(){
pcout << "Assemble" << endl;
// Set scope of the timer
TimerOutput::Scope t(timer, "Assembling");
system_matrix = 0;
system_rhs = 0;
// Define quadrature rule and degree
QGauss<DIM-1> face_quadrature_formula(params::degree+2);
FEFaceValues<DIM> fe_face_values (fe, face_quadrature_formula,
update_values |
update_normal_vectors |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_face_q_points = face_quadrature_formula.size();
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);
typename DoFHandler<DIM>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell) {
// If not in this processor, move on
if(!cell->is_locally_owned()) continue;
cell_matrix = 0;
cell_rhs = 0;
for (auto face_num=0u; face_num<GeometryInfo<DIM>::faces_per_cell; ++face_num){
if (cell->face(face_num)->at_boundary()
&& (cell->face(face_num)->boundary_id() == 1 ) ){
// Update face values
fe_face_values.reinit (cell, face_num);
for (auto q=0u; q < n_face_q_points; ++q){
for (auto i=0u; i<dofs_per_cell; ++i){
cell_rhs(i) += fe_face_values.shape_value(i, q) *
2 *
fe_face_values.JxW(q);
}
}
}// end if at boundary
}// end face
// local-to-global
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global(cell_matrix, cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
// end local-to-global
}
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
}
void spGeom::output_results(){
pcout << "Output" << endl;
TimerOutput::Scope t(timer, "Output results");
// Set solution names
std::vector<string> solution_names (DIM, "displacements");
solution_names.push_back ("pressure");
// Set the type of value DataOut class should expect
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(DIM,
DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation.push_back (DataComponentInterpretation::component_is_scalar);
DataOut<DIM> data_out;
data_out.attach_dof_handler (dof_handler);
LA::MPI::BlockVector tmp;
tmp.reinit (owned_partitioning,
relevant_partitioning,
mpi_communicator);
tmp = system_rhs;
data_out.add_data_vector (tmp,
solution_names,
DataOut<DIM>::type_dof_data,
data_component_interpretation);
data_out.build_patches ();
const string filename = ("rhs-" +
Utilities::int_to_string
(triangulation.locally_owned_subdomain(), 4) + ".vtu");
std::ofstream output (filename.c_str());
data_out.write_vtu (output);
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
std::vector<string> filenames;
for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
filenames.push_back (string("rhs-") +
Utilities::int_to_string(i, 4) + ".vtu");
const string pvtu_master_filename = ("rhs.pvtu");
std::ofstream pvtu_master (pvtu_master_filename.c_str());
data_out.write_pvtu_record (pvtu_master, filenames);
}
}
void spGeom::run(){
makeGrid();
setupDofs();
assemble();
output_results ();
}
int main (int argc, char** argv)
{
using namespace dealii;
using namespace std;
try{
// Initialize MPI
Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
string filename = "gia.log";
ofstream logFile(filename.c_str());
// Attach deallog to process output
deallog.attach(logFile);
deallog.depth_console (0);
spGeom ex;
ex.run ();
}
catch (exception &exc){
cerr << setw(30) << setfill('-') << "" << endl;
cerr << "Exception on processing: " << endl
<< exc.what() << endl
<< "Aborting!" << endl
<< setw(30) << setfill('-') << "" << endl;
return EXIT_FAILURE;
}
catch (...){
cerr << setw(30) << setfill('-') << "" << endl;
cerr << "Unknown exception!" << endl
<< "Aborting!" << endl
<< setw(30) << setfill('-') << "" << endl;
return EXIT_FAILURE;
}
return 0;
}
