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

Reply via email to