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; }