From 44c9c47ded3cdcf682268c09abe1dc7b1daa4099 Mon Sep 17 00:00:00 2001 From: Doug Shi-Dong Date: Fri, 4 Oct 2019 22:22:47 -0400 Subject: [PATCH] Implementat HighOrderGrid with all existing test Had to slightly change the Gaussian bump test to ensure convergence by increasing the initial grid size. As a result, we do one less grid level to keep it tractable on a personal computer. Also fixed a bug regarding the DoFCellAccessor being used with the Lagrange FiniteElement. --- src/dg/dg.cpp | 314 +++++------------- src/dg/dg.h | 51 ++- src/dg/high_order_grid.cpp | 42 ++- src/dg/high_order_grid.h | 36 +- src/dg/strong_dg.cpp | 2 - src/testing/euler_cylinder.cpp | 4 +- src/testing/euler_gaussian_bump.cpp | 18 +- src/testing/grid_study.cpp | 4 +- .../2d_euler_gaussian_bump.prm | 4 +- .../unit_tests/grid/high_order_grid_test.cpp | 7 +- 10 files changed, 192 insertions(+), 290 deletions(-) diff --git a/src/dg/dg.cpp b/src/dg/dg.cpp index 3e24192c6..b9099ec87 100644 --- a/src/dg/dg.cpp +++ b/src/dg/dg.cpp @@ -9,6 +9,8 @@ #include #include +#include + #include #include #include @@ -46,136 +48,6 @@ namespace PHiLiP { template using Triangulation = dealii::parallel::distributed::Triangulation; #endif -//template -// > -//void -//get_position_vector2(const dealii::hp::DoFHandler &dh, VectorType &vector) -//{ -// AssertDimension(vector.size(), dh.n_dofs()); -// const dealii::hp::FECollection &fe_coll = dh.get_fe_collection(); -// -// dealii::hp::MappingCollection mapping_coll; -// dealii::hp::QCollection quad_coll; -// for (int i=0; i finite_element = fe_coll[i]; -// -// Assert(finite_element.is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); -// assert(finite_element.has_support_points()); -// AssertDimension(finite_element.dofs_per_cell, finite_element.get_unit_support_points().size()); -// -// const dealii::Quadrature quadrature(finite_element.get_unit_support_points()); -// const dealii::MappingQ mapping(finite_element.degree); -// -// mapping_coll.push_back(mapping); -// quad_coll.push_back(quadrature); -// } -// -// dealii::hp::FEValues fe_values_coll (mapping_coll, fe_coll, quadrature_coll, update_quadrature_points); -// -// // Construct default fe_mask; -// const ComponentMask fe_mask(//mask.size() ? mask : -// ComponentMask(fe_coll.get_nonzero_components(0).size(), true)); -// -// AssertDimension(fe_mask.size(), fe_coll.get_nonzero_components(0).size()); -// -// std::vector fe_to_real(fe_mask.size(), -// numbers::invalid_unsigned_int); -// unsigned int size = 0; -// for (unsigned int i = 0; i < fe_mask.size(); ++i) { -// if (fe_mask[i]) fe_to_real[i] = size++; -// } -// Assert(size == spacedim, -// ExcMessage("The Component Mask you provided is invalid. It has to select exactly spacedim entries.")); -// -// std::vector dof_indices(fe_coll[fe_coll.size()-1].dofs_per_cell); -// for (const auto &cell : dh.active_cell_iterators()) { -// if (cell->is_locally_owned()) { -// const unsigned int fe_index = cell->active_fe_index(); -// fe_values_coll.reinit (cell, fe_index, fe_index, fe_index); -// dof_indices.resize(fe_coll[fe_index].dofs_per_cell); -// cell->get_dof_indices(dof_indices); -// const std::vector> &points = fe_v.get_quadrature_points(); -// for (unsigned int q = 0; q < points.size(); ++q) { -// -// const unsigned int comp = fe_coll.system_to_component_index(q).first; -// if (fe_mask[comp]) { -// dealii::internal::ElementAccess:: -// set(points[q][fe_to_real[comp]], dof_indices[q], vector); -// } -// } -// } -// } -//} - -//template -// > -//void -//get_position_vector2(const dealii::hp::DoFHandler &dof_handler, VectorType &vector) -//{ -// AssertDimension(vector.size(), dof_handler.n_dofs()); -// const dealii::hp::FECollection &fe_coll = dof_handler.get_fe_collection(); -// -// dealii::hp::MappingCollection mapping_coll; -// dealii::hp::QCollection quad_coll; -// for (unsigned int i=0; i &finite_element = fe_coll[i]; -// -// Assert(finite_element.is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); -// assert(finite_element.has_support_points()); -// AssertDimension(finite_element.dofs_per_cell, finite_element.get_unit_support_points().size()); -// -// const dealii::Quadrature quadrature(finite_element.get_unit_support_points()); -// const dealii::MappingQ mapping(finite_element.degree); -// -// mapping_coll.push_back(mapping); -// quad_coll.push_back(quadrature); -// } -// -// dealii::hp::FEValues fe_values_coll (mapping_coll, fe_coll, quad_coll, dealii::update_quadrature_points); -// -// // Construct default fe_mask; -// const dealii::ComponentMask fe_mask(//mask.size() ? mask : -// dealii::ComponentMask(fe_coll[0].get_nonzero_components(0).size(), true)); -// -// AssertDimension(fe_mask.size(), fe_coll[0].get_nonzero_components(0).size()); -// -// std::vector fe_to_real(fe_mask.size(), -// dealii::numbers::invalid_unsigned_int); -// unsigned int size = 0; -// for (unsigned int i = 0; i < fe_mask.size(); ++i) { -// if (fe_mask[i]) fe_to_real[i] = size++; -// } -// Assert(size == spacedim, -// ExcMessage("The Component Mask you provided is invalid. It has to select exactly spacedim entries.")); -// -// std::vector dof_indices(fe_coll[fe_coll.size()-1].dofs_per_cell); -// for (const auto &cell : dof_handler.active_cell_iterators()) { -// if (cell->is_locally_owned()) { -// const unsigned int fe_index = cell->active_fe_index(); -// fe_values_coll.reinit (cell, fe_index, fe_index, fe_index); -// const dealii::FEValues &fe_values_volume = fe_values_coll.get_present_fe_values(); -// -// dof_indices.resize(fe_coll[fe_index].dofs_per_cell); -// cell->get_dof_indices(dof_indices); -// -// const std::vector> &points = fe_values_volume.get_quadrature_points(); -// for (unsigned int q = 0; q < points.size(); ++q) { -// -// const unsigned int comp = fe_values_volume.get_fe().system_to_component_index(q).first; -// if (fe_mask[comp]) { -// dealii::internal::ElementAccess:: -// set(points[q][fe_to_real[comp]], dof_indices[q], vector); -// } -// } -// } -// } -//} - - // DGFactory *********************************************************************** template std::shared_ptr< DGBase > @@ -223,7 +95,7 @@ ::create_discontinuous_galerkin( // DGBase *************************************************************************** template -DGBase::DGBase( // @suppress("Class members should be properly initialized") +DGBase::DGBase( const int nstate_input, const Parameters::AllParameters *const parameters_input, const unsigned int max_degree_input, @@ -232,7 +104,7 @@ DGBase::DGBase( // @suppress("Class members should be properly initial { } template -DGBase::DGBase( // @suppress("Class members should be properly initialized") +DGBase::DGBase( const int nstate_input, const Parameters::AllParameters *const parameters_input, const unsigned int max_degree_input, @@ -242,26 +114,18 @@ DGBase::DGBase( // @suppress("Class members should be properly initial , nstate(nstate_input) , max_degree(max_degree_input) , triangulation(triangulation_input) - //, mapping_collection(std::get<0>(collection_tuple)) , fe_collection(std::get<0>(collection_tuple)) - , fe_grid(std::get<1>(collection_tuple),dim) - , volume_quadrature_collection(std::get<2>(collection_tuple)) - , face_quadrature_collection(std::get<3>(collection_tuple)) - , oned_quadrature_collection(std::get<4>(collection_tuple)) - , fe_collection_lagrange(std::get<5>(collection_tuple)) + , volume_quadrature_collection(std::get<1>(collection_tuple)) + , face_quadrature_collection(std::get<2>(collection_tuple)) + , oned_quadrature_collection(std::get<3>(collection_tuple)) + , fe_collection_lagrange(std::get<4>(collection_tuple)) , dof_handler(*triangulation) - , dof_handler_grid(*triangulation) - , high_order_grid_mapping(dof_handler_grid, high_order_grid) + , high_order_grid(all_parameters, max_degree_input+1, triangulation) , mpi_communicator(MPI_COMM_WORLD) , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) { dof_handler.initialize(*triangulation, fe_collection); - dof_handler_grid.initialize(*triangulation, fe_grid); - int minimum_degree = (all_parameters->use_collocated_nodes==true) ? 1 : 0; - for (unsigned int degree=minimum_degree; degree<=max_degree; ++degree) { - mapping_collection.push_back(high_order_grid_mapping); // Can't just return FESystem. For some reason the copy constructor is deleted - } set_all_cells_fe_degree (fe_collection.size()-1); // Always sets it to the maximum degree @@ -271,8 +135,6 @@ template std::tuple< //dealii::hp::MappingCollection, // Mapping dealii::hp::FECollection, // Solution FE - //dealii::hp::FECollection, // Grid FE - dealii::FE_Q, // Grid FE dealii::hp::QCollection, // Volume quadrature dealii::hp::QCollection, // Face quadrature dealii::hp::QCollection<1>, // 1D quadrature for strong form @@ -281,8 +143,6 @@ DGBase::create_collection_tuple(const unsigned int max_degree, const i { //dealii::hp::MappingCollection mapping_coll; dealii::hp::FECollection fe_coll; - //dealii::hp::FECollection fe_coll_grid; - const dealii::FE_Q feq_grid(max_degree+1); // The grid must be at least p1. A p0 solution required a p1 grid. dealii::hp::QCollection volume_quad_coll; dealii::hp::QCollection face_quad_coll; dealii::hp::QCollection<1> oned_quad_coll; @@ -294,53 +154,40 @@ DGBase::create_collection_tuple(const unsigned int max_degree, const i //const dealii::MappingQ mapping(degree+1, true); //const dealii::MappingManifold mapping; //mapping_coll.push_back(mapping); - //mapping_coll.push_back(high_order_grid_mapping); // Can't just return FESystem. For some reason the copy constructor is deleted // Solution FECollection const dealii::FE_DGQ fe_dg(degree); const dealii::FESystem fe_system(fe_dg, nstate); fe_coll.push_back (fe_system); - //// Grid FECollection uses Lagrange polynomials - //const dealii::FE_Q fe_q(degree+1); // The grid must be at least p1. A p0 solution required a p1 grid. - //const dealii::FESystem fe_system_q(fe_q, dim); - //fe_coll_grid.push_back (fe_system_q); - dealii::Quadrature<1> oned_quad(degree+1); dealii::Quadrature volume_quad(degree+1); dealii::Quadrature face_quad(degree+1); //removed const - if (parameters_input->use_collocated_nodes) - { - dealii::QGaussLobatto<1> oned_quad_Gauss_Lobatto (degree+1); - dealii::QGaussLobatto vol_quad_Gauss_Lobatto (degree+1); - oned_quad = oned_quad_Gauss_Lobatto; - volume_quad = vol_quad_Gauss_Lobatto; - - if(dim == 1) - { - dealii::QGauss face_quad_Gauss_Legendre (degree+1); - face_quad = face_quad_Gauss_Legendre; - } - else - { - dealii::QGaussLobatto face_quad_Gauss_Lobatto (degree+1); - face_quad = face_quad_Gauss_Lobatto; - } - + if (parameters_input->use_collocated_nodes) { + dealii::QGaussLobatto<1> oned_quad_Gauss_Lobatto (degree+1); + dealii::QGaussLobatto vol_quad_Gauss_Lobatto (degree+1); + oned_quad = oned_quad_Gauss_Lobatto; + volume_quad = vol_quad_Gauss_Lobatto; - } - else + if(dim == 1) { - dealii::QGauss<1> oned_quad_Gauss_Legendre (degree+1); - dealii::QGauss vol_quad_Gauss_Legendre (degree+1); dealii::QGauss face_quad_Gauss_Legendre (degree+1); - oned_quad = oned_quad_Gauss_Legendre; - volume_quad = vol_quad_Gauss_Legendre; face_quad = face_quad_Gauss_Legendre; } - // - + else + { + dealii::QGaussLobatto face_quad_Gauss_Lobatto (degree+1); + face_quad = face_quad_Gauss_Lobatto; + } + } else { + dealii::QGauss<1> oned_quad_Gauss_Legendre (degree+1); + dealii::QGauss vol_quad_Gauss_Legendre (degree+1); + dealii::QGauss face_quad_Gauss_Legendre (degree+1); + oned_quad = oned_quad_Gauss_Legendre; + volume_quad = vol_quad_Gauss_Legendre; + face_quad = face_quad_Gauss_Legendre; + } volume_quad_coll.push_back (volume_quad); face_quad_coll.push_back (face_quad); @@ -349,7 +196,7 @@ DGBase::create_collection_tuple(const unsigned int max_degree, const i dealii::FE_DGQArbitraryNodes lagrange_poly(oned_quad); fe_coll_lagr.push_back (lagrange_poly); } - return std::make_tuple(fe_coll, feq_grid, volume_quad_coll, face_quad_coll, oned_quad_coll, fe_coll_lagr); + return std::make_tuple(fe_coll, volume_quad_coll, face_quad_coll, oned_quad_coll, fe_coll_lagr); } @@ -362,10 +209,6 @@ void DGBase::set_all_cells_fe_degree ( const unsigned int degree ) if (cell->is_locally_owned()) cell->set_future_fe_index (degree); } - //for (auto cell = dof_handler_grid.begin_active(); cell != dof_handler_grid.end(); ++cell) - //{ - // if (cell->is_locally_owned()) cell->set_future_fe_index (degree); - //} triangulation->execute_coarsening_and_refinement(); } @@ -376,19 +219,16 @@ template DGBase::~DGBase () { dof_handler.clear (); - dof_handler_grid.clear (); } // template // void DGBase::set_triangulation(Triangulation *triangulation_input) // { // // dof_handler.clear(); -// // dof_handler_grid.clear(); // // // triangulation = triangulation_input; // // // dof_handler.initialize(*triangulation, fe_collection); -// // dof_handler_grid.initialize(*triangulation, fe_grid); // // //set_all_cells_fe_degree (fe_collection.size()-1); // Always sets it to the maximum degree // @@ -408,6 +248,12 @@ void DGBase::assemble_residual (const bool compute_dRdW) std::vector current_dofs_indices(max_dofs_per_cell); std::vector neighbor_dofs_indices(max_dofs_per_cell); + //dealii::hp::MappingCollection mapping_collection(*(high_order_grid.mapping_fe_field)); + //const dealii::MappingManifold mapping; + //const dealii::MappingQ mapping(max_degree+1); + const auto mapping = (*(high_order_grid.mapping_fe_field)); + dealii::hp::MappingCollection mapping_collection(mapping); + dealii::hp::FEValues fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags); ///< FEValues of volume. dealii::hp::FEFaceValues fe_values_collection_face_int (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of interior face. dealii::hp::FEFaceValues fe_values_collection_face_ext (mapping_collection, fe_collection, face_quadrature_collection, this->neighbor_face_update_flags); ///< FEValues of exterior face. @@ -424,7 +270,9 @@ void DGBase::assemble_residual (const bool compute_dRdW) n_cell_visited++; // Current reference element related to this physical cell + const unsigned int mapping_index = 0; const unsigned int fe_index_curr_cell = current_cell->active_fe_index(); + const unsigned int quad_index = fe_index_curr_cell; const dealii::FESystem ¤t_fe_ref = fe_collection[fe_index_curr_cell]; const unsigned int curr_cell_degree = current_fe_ref.tensor_degree(); const unsigned int n_dofs_curr_cell = current_fe_ref.n_dofs_per_cell(); @@ -437,10 +285,13 @@ void DGBase::assemble_residual (const bool compute_dRdW) current_cell->get_dof_indices (current_dofs_indices); // fe_values_collection.reinit(current_cell, quad_collection_index, mapping_collection_index, fe_collection_index) - fe_values_collection_volume.reinit (current_cell, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_volume.reinit (current_cell, quad_index, mapping_index, fe_index_curr_cell); const dealii::FEValues &fe_values_volume = fe_values_collection_volume.get_present_fe_values(); - if (!(all_parameters->use_weak_form)) fe_values_collection_volume_lagrange.reinit (current_cell, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + + dealii::TriaIterator> cell_iterator = static_cast> > (current_cell); + //if (!(all_parameters->use_weak_form)) fe_values_collection_volume_lagrange.reinit (current_cell, quad_index, mapping_index, fe_index_curr_cell); + fe_values_collection_volume_lagrange.reinit (cell_iterator, quad_index, mapping_index, fe_index_curr_cell); const dealii::FEValues &fe_values_lagrange = fe_values_collection_volume_lagrange.get_present_fe_values(); if ( compute_dRdW ) { assemble_volume_terms_implicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); @@ -461,7 +312,7 @@ void DGBase::assemble_residual (const bool compute_dRdW) n_face_visited++; - fe_values_collection_face_int.reinit (current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); if(current_face->at_boundary() && all_parameters->use_periodic_bc == true && dim == 1) //using periodic BCs (for 1d) { @@ -469,7 +320,7 @@ void DGBase::assemble_residual (const bool compute_dRdW) //int cell_index = current_cell->index(); if (cell_index == 0 && iface == 0) { - fe_values_collection_face_int.reinit(current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_face_int.reinit(current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); neighbor_cell = dof_handler.begin_active(); for (unsigned int i = 0 ; i < triangulation->n_active_cells() - 1; ++i) { @@ -477,24 +328,27 @@ void DGBase::assemble_residual (const bool compute_dRdW) } neighbor_cell->get_dof_indices(neighbor_dofs_indices); const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); + const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; + const unsigned int mapping_index_neigh_cell = 0; - fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1,fe_index_neigh_cell,fe_index_neigh_cell,fe_index_neigh_cell); + fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1,quad_index_neigh_cell,mapping_index_neigh_cell,fe_index_neigh_cell); } else if (cell_index == (int) triangulation->n_active_cells() - 1 && iface == 1) { - fe_values_collection_face_int.reinit(current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_face_int.reinit(current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); neighbor_cell = dof_handler.begin_active(); neighbor_cell->get_dof_indices(neighbor_dofs_indices); - const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); - fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1, fe_index_neigh_cell, fe_index_neigh_cell, fe_index_neigh_cell); //not sure how changing the face number would work in dim!=1-dimensions. + const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); + const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; + const unsigned int mapping_index_neigh_cell = 0; + fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); //not sure how changing the face number would work in dim!=1-dimensions. } //std::cout << "cell " << current_cell->index() << "'s " << iface << "th face has neighbour: " << neighbor_cell->index() << std::endl; const int neighbor_face_no = (iface ==1) ? 0:1; const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); - const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); @@ -575,6 +429,8 @@ void DGBase::assemble_residual (const bool compute_dRdW) // Get information about neighbor cell const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); + const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; + const unsigned int mapping_index_neigh_cell = 0; const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); @@ -586,9 +442,9 @@ void DGBase::assemble_residual (const bool compute_dRdW) neighbor_dofs_indices.resize(n_dofs_neigh_cell); neighbor_cell->get_dof_indices (neighbor_dofs_indices); - fe_values_collection_face_int.reinit (current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, fe_index_neigh_cell, fe_index_neigh_cell, fe_index_neigh_cell); + fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; @@ -652,6 +508,8 @@ void DGBase::assemble_residual (const bool compute_dRdW) // Get information about neighbor cell const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); + const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; + const unsigned int mapping_index_neigh_cell = 0; const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); @@ -662,10 +520,10 @@ void DGBase::assemble_residual (const bool compute_dRdW) neighbor_dofs_indices.resize(n_dofs_neigh_cell); neighbor_cell->get_dof_indices (neighbor_dofs_indices); - fe_values_collection_subface.reinit (current_cell, iface, subface_no, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_subface.reinit (current_cell, iface, subface_no, quad_index, mapping_index, fe_index_curr_cell); const dealii::FESubfaceValues &fe_values_face_int = fe_values_collection_subface.get_present_fe_values(); - fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, fe_index_neigh_cell, fe_index_neigh_cell, fe_index_neigh_cell); + fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; @@ -737,6 +595,8 @@ void DGBase::assemble_residual (const bool compute_dRdW) // Get information about neighbor cell const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); + const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; + const unsigned int mapping_index_neigh_cell = 0; const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); @@ -748,9 +608,9 @@ void DGBase::assemble_residual (const bool compute_dRdW) neighbor_dofs_indices.resize(n_dofs_neigh_cell); neighbor_cell->get_dof_indices (neighbor_dofs_indices); - fe_values_collection_face_int.reinit (current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, fe_index_neigh_cell, fe_index_neigh_cell, fe_index_neigh_cell); + fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; @@ -868,7 +728,9 @@ void DGBase::output_results_vtk (const unsigned int cycle)// const const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); - data_out.build_patches (mapping_collection[mapping_collection.size()-1]); + //data_out.build_patches (mapping_collection[mapping_collection.size()-1]); + data_out.build_patches(*(high_order_grid.mapping_fe_field), max_degree, dealii::DataOut>::CurvedCellRegion::curved_inner_cells); + //data_out.build_patches(*(high_order_grid.mapping_fe_field), fe_collection.size(), dealii::DataOut::CurvedCellRegion::curved_inner_cells); std::string filename = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D-"; filename += dealii::Utilities::int_to_string(cycle, 4) + "."; filename += dealii::Utilities::int_to_string(iproc, 4); @@ -901,28 +763,20 @@ void DGBase::allocate_system () // system matrices and vectors. dof_handler.distribute_dofs(fe_collection); - dof_handler_grid.distribute_dofs(fe_grid); - - // High order grid initialization - locally_owned_dofs_grid = dof_handler_grid.locally_owned_dofs(); - dealii::DoFTools::extract_locally_relevant_dofs(dof_handler_grid, ghost_dofs_grid); - locally_relevant_dofs_grid = ghost_dofs_grid; - ghost_dofs_grid.subtract_set(locally_owned_dofs_grid); - high_order_grid.reinit(locally_owned_dofs_grid, ghost_dofs_grid, mpi_communicator); - - //dealii::VectorTools::get_position_vector2(dof_handler_grid, high_order_grid); - //get_position_vector2(dof_handler_grid, high_order_grid); - dealii::MappingFEField, dealii::DoFHandler > map(dof_handler_grid, high_order_grid); + + //dealii::MappingFEField, dealii::DoFHandler> mapping = high_order_grid.get_MappingFEField(); + //dealii::MappingFEField, dealii::DoFHandler> mapping = *(high_order_grid.mapping_fe_field); + //int minimum_degree = (all_parameters->use_collocated_nodes==true) ? 1 : 0; - ////int current_fe_index = 0; + //int current_fe_index = 0; //for (unsigned int degree=minimum_degree; degree<=max_degree; ++degree) { - // mapping_collection.push_back(map); - // //if(current_fe_index <= mapping_collection.size()) { - // // mapping_collection.push_back(map); - // //} else { - // // mapping_collection[current_fe_index] = map; - // //} - // //current_fe_index++; + // //mapping_collection.push_back(mapping); + // if(current_fe_index <= mapping_collection.size()) { + // mapping_collection.push_back(mapping); + // } else { + // mapping_collection[current_fe_index] = std::shared_ptr>(mapping.clone()); + // } + // current_fe_index++; //} // Solution and RHS @@ -994,12 +848,20 @@ void DGBase::evaluate_mass_matrices (bool do_inverse_mass_matrix) //} //pcout << "AFter reinit" << std::endl; + //dealii::hp::MappingCollection mapping_collection(*(high_order_grid.mapping_fe_field)); + //const dealii::MappingManifold mapping; + //const dealii::MappingQ mapping(max_degree+1); + const auto mapping = (*(high_order_grid.mapping_fe_field)); + dealii::hp::MappingCollection mapping_collection(mapping); + dealii::hp::FEValues fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags); ///< FEValues of volume. for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) { if (!cell->is_locally_owned()) continue; + const unsigned int mapping_index = 0; const unsigned int fe_index_curr_cell = cell->active_fe_index(); + const unsigned int quad_index = fe_index_curr_cell; // Current reference element related to this physical cell const dealii::FESystem ¤t_fe_ref = fe_collection[fe_index_curr_cell]; @@ -1008,7 +870,7 @@ void DGBase::evaluate_mass_matrices (bool do_inverse_mass_matrix) dealii::FullMatrix local_mass_matrix(n_dofs_cell); - fe_values_collection_volume.reinit (cell, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell); + fe_values_collection_volume.reinit (cell, quad_index, mapping_index, fe_index_curr_cell); const dealii::FEValues &fe_values_volume = fe_values_collection_volume.get_present_fe_values(); for (unsigned int itest=0; itest +#include "dg/high_order_grid.h" #include "physics/physics.h" #include "numerical_flux/numerical_flux.h" #include "parameters/all_parameters.h" @@ -100,9 +101,6 @@ class DGBase using MassiveCollectionTuple = std::tuple< //dealii::hp::MappingCollection, // Mapping dealii::hp::FECollection, // Solution FE - //dealii::hp::FECollection, // Grid FE - //dealii::FESystem, // Grid FE // Can't just return FESystem. For some reason the copy constructor is deleted - dealii::FE_Q, // Grid FE dealii::hp::QCollection, // Volume quadrature dealii::hp::QCollection, // Face quadrature dealii::hp::QCollection<1>, // 1D quadrature for strong form @@ -129,7 +127,7 @@ class DGBase /** Since we are interested in performing mesh movement for optimization purposes, * this is not a constant member variables. */ - dealii::hp::MappingCollection mapping_collection; + //dealii::hp::MappingCollection mapping_collection; void set_all_cells_fe_degree ( const unsigned int degree ); @@ -219,9 +217,6 @@ class DGBase */ dealii::LinearAlgebra::distributed::Vector solution; - /// Current nodal coefficients of the high-order grid. - dealii::LinearAlgebra::distributed::Vector high_order_grid; - void initialize_manufactured_solution (); ///< Virtual function defined in DG void output_results_vtk (const unsigned int ith_grid); ///< Output solution @@ -271,13 +266,30 @@ class DGBase * Therefore, every grid/cell will use the maximal polynomial mapping regardless of the solution order. */ //const dealii::hp::FECollection fe_collection_grid; - const dealii::FESystem fe_grid; + //const dealii::FESystem fe_grid; dealii::hp::QCollection volume_quadrature_collection; dealii::hp::QCollection face_quadrature_collection; dealii::hp::QCollection<1> oned_quadrature_collection; +protected: + /// Lagrange basis used in strong form + /** This is a collection of scalar Lagrange bases */ + const dealii::hp::FECollection fe_collection_lagrange; +public: + /// Degrees of freedom handler + /* Allows us to iterate over the finite elements' degrees of freedom. + * Note that since we are not using FESystem, we need to multiply + * the index by a factor of "nstate" + * + * Must be defined after fe_dg since it is a subscriptor of fe_dg. + * Destructor are called in reverse order in which they appear in class definition. + */ + dealii::hp::DoFHandler dof_handler; + + /// High order grid that will provide the MappingFEField + HighOrderGrid high_order_grid; protected: @@ -370,31 +382,8 @@ class DGBase const dealii::UpdateFlags neighbor_face_update_flags = dealii::update_values | dealii::update_gradients; - /// Lagrange basis used in strong form - /** This is a collection of scalar Lagrange bases */ - const dealii::hp::FECollection fe_collection_lagrange; - //dealii::hp::FEValues fe_values_collection_volume_lagrange; - - -public: - /// Degrees of freedom handler - /* Allows us to iterate over the finite elements' degrees of freedom. - * Note that since we are not using FESystem, we need to multiply - * the index by a factor of "nstate" - * - * Must be defined after fe_dg since it is a subscriptor of fe_dg. - * Destructor are called in reverse order in which they appear in class definition. - */ - dealii::hp::DoFHandler dof_handler; - - /// Degrees of freedom handler for the high-order grid - dealii::DoFHandler dof_handler_grid; protected: - /// Main underlying mapping responsible for all mapping evaluations - //dealii::MappingFEField, dealii::hp::DoFHandler > high_order_grid_mapping; - dealii::MappingFEField, dealii::DoFHandler > high_order_grid_mapping; - MPI_Comm mpi_communicator; ///< MPI communicator dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 private: diff --git a/src/dg/high_order_grid.cpp b/src/dg/high_order_grid.cpp index 23ce65b87..500a37b47 100644 --- a/src/dg/high_order_grid.cpp +++ b/src/dg/high_order_grid.cpp @@ -23,6 +23,10 @@ HighOrderGrid::HighOrderGrid( , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) { allocate(); + const dealii::ComponentMask mask(dim, true); + dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask); + nodes.update_ghost_values(); + mapping_fe_field = std::make_shared< dealii::MappingFEField > (dof_handler_grid,nodes,mask); } template @@ -32,23 +36,28 @@ HighOrderGrid::allocate() dof_handler_grid.initialize(*triangulation, fe_system); dof_handler_grid.distribute_dofs(fe_system); +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + nodes.reinit(dof_handler_grid.n_dofs()); +#else locally_owned_dofs_grid = dof_handler_grid.locally_owned_dofs(); dealii::DoFTools::extract_locally_relevant_dofs(dof_handler_grid, ghost_dofs_grid); locally_relevant_dofs_grid = ghost_dofs_grid; ghost_dofs_grid.subtract_set(locally_owned_dofs_grid); nodes.reinit(locally_owned_dofs_grid, ghost_dofs_grid, mpi_communicator); +#endif } -template -dealii::MappingFEField -HighOrderGrid::get_MappingFEField() { - const dealii::ComponentMask mask(dim, true); - dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask); - - dealii::MappingFEField mapping(dof_handler_grid,nodes,mask); - - return mapping; -} +//template +//dealii::MappingFEField +//HighOrderGrid::get_MappingFEField() { +// const dealii::ComponentMask mask(dim, true); +// dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask); +// +// dealii::MappingFEField mapping(dof_handler_grid,nodes,mask); +// +// return mapping; +//} +// template @@ -56,17 +65,30 @@ void HighOrderGrid::prepare_for_coarsening_a old_nodes = nodes; old_nodes.update_ghost_values(); +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + solution_transfer.prepare_for_coarsening_and_refinement(old_nodes); +#else solution_transfer.prepare_for_coarsening_and_refinement(old_nodes); +#endif } template void HighOrderGrid::execute_coarsening_and_refinement() { allocate(); +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + solution_transfer.interpolate(old_nodes, nodes); +#else solution_transfer.interpolate(nodes); +#endif nodes.update_ghost_values(); + mapping_fe_field = std::make_shared< dealii::MappingFEField > (dof_handler_grid,nodes); } //template class HighOrderGrid; +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D +template class HighOrderGrid, dealii::DoFHandler>; +#else template class HighOrderGrid, dealii::DoFHandler>; +#endif } // namespace PHiLiP diff --git a/src/dg/high_order_grid.h b/src/dg/high_order_grid.h index 75bd885d6..30ea31638 100644 --- a/src/dg/high_order_grid.h +++ b/src/dg/high_order_grid.h @@ -8,6 +8,7 @@ #include +#include #include #include @@ -22,10 +23,24 @@ namespace PHiLiP { * nodes vector. The dof_handler_grid is used to access and loop through those nodes. * This will especially be useful when performing shape optimization, where the surface and volume * nodes will need to be displaced. + * Note that there are a lot of pre-processor statements, and that is because the SolutionTransfer class + * and the Vector class act quite differently between serial and parallel implementation. Hopefully, + * deal.II will change this one day such that we have one interface for both. */ +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D +template , typename DoFHandlerType = dealii::DoFHandler> +#else template , typename DoFHandlerType = dealii::DoFHandler> +#endif class HighOrderGrid { +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + using Vector = dealii::Vector; + using SolutionTransfer = dealii::SolutionTransfer, dealii::DoFHandler>; +#else + using Vector = dealii::LinearAlgebra::distributed::Vector; + using SolutionTransfer = dealii::parallel::distributed::SolutionTransfer, dealii::DoFHandler>; +#endif public: /// Principal constructor that will call delegated constructor. HighOrderGrid(const Parameters::AllParameters *const parameters_input, const unsigned int max_degree, dealii::Triangulation *const triangulation_input); @@ -36,6 +51,9 @@ class HighOrderGrid /// Return a MappingFEField that corresponds to the current node locations dealii::MappingFEField get_MappingFEField(); + /// Return a MappingFEField that corresponds to the current node locations + void update_MappingFEField(); + const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters /// Maximum degree of the geometry polynomial representing the grid. @@ -47,13 +65,13 @@ class HighOrderGrid dealii::DoFHandler dof_handler_grid; /// Current nodal coefficients of the high-order grid. - dealii::LinearAlgebra::distributed::Vector nodes; + Vector nodes; /// List of surface points - dealii::LinearAlgebra::distributed::Vector surface_nodes; + Vector surface_nodes; /// RBF mesh deformation - To be done - void deform_mesh(dealii::LinearAlgebra::distributed::Vector surface_displacements); + void deform_mesh(Vector surface_displacements); /// Evaluate cell metric Jacobian /** The metric Jacobian is given by the gradient of the physical location @@ -84,18 +102,26 @@ class HighOrderGrid /// Using system of polynomials to represent the x, y, and z directions. const dealii::FESystem fe_system; + + /** MappingFEField that will provide the polynomial-based grid. + * It is a shared smart pointer because the constructor requires the dof_handler_grid to be properly initialized. + * See discussion in the following thread: + * https://stackoverflow.com/questions/7557153/defining-an-object-without-calling-its-constructor-in-c + */ + std::shared_ptr> mapping_fe_field; + dealii::IndexSet locally_owned_dofs_grid; ///< Locally own degrees of freedom for the grid dealii::IndexSet ghost_dofs_grid; ///< Locally relevant ghost degrees of freedom for the grid dealii::IndexSet locally_relevant_dofs_grid; ///< Union of locally owned degrees of freedom and relevant ghost degrees of freedom for the grid protected: /// Used for the SolutionTransfer when performing grid adaptation. - dealii::LinearAlgebra::distributed::Vector old_nodes; + Vector old_nodes; /** Transfers the coarse curved curve onto the fine curved grid. * Used in prepare_for_coarsening_and_refinement() and execute_coarsening_and_refinement() */ - dealii::parallel::distributed::SolutionTransfer, dealii::DoFHandler> solution_transfer; + SolutionTransfer solution_transfer; MPI_Comm mpi_communicator; ///< MPI communicator dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 diff --git a/src/dg/strong_dg.cpp b/src/dg/strong_dg.cpp index a3ae702fc..9bfb24c87 100644 --- a/src/dg/strong_dg.cpp +++ b/src/dg/strong_dg.cpp @@ -581,8 +581,6 @@ void DGStrong::assemble_volume_terms_explicit( // Convective // Now minus such 2 integrations by parts - assert(JxW[iquad] - fe_values_lagrange.JxW(iquad) < 1e-14); - rhs = rhs - fe_values_vol.shape_value_component(itest,iquad,istate) * flux_divergence[iquad][istate] * JxW[iquad]; //// Diffusive diff --git a/src/testing/euler_cylinder.cpp b/src/testing/euler_cylinder.cpp index be96398ec..46952071e 100644 --- a/src/testing/euler_cylinder.cpp +++ b/src/testing/euler_cylinder.cpp @@ -206,7 +206,9 @@ ::run_test () const dealii::LinearAlgebra::distributed::Vector old_solution(dg->solution); dealii::parallel::distributed::SolutionTransfer, dealii::hp::DoFHandler> solution_transfer(dg->dof_handler); solution_transfer.prepare_for_coarsening_and_refinement(old_solution); + dg->high_order_grid.prepare_for_coarsening_and_refinement(); grid.refine_global (1); + dg->high_order_grid.execute_coarsening_and_refinement(); dg->allocate_system (); dg->solution.zero_out_ghosts(); solution_transfer.interpolate(dg->solution); @@ -236,7 +238,7 @@ ::run_test () const // Overintegrate the error to make sure there is not integration error in the error estimate int overintegrate = 10; dealii::QGauss quad_extra(dg->max_degree+1+overintegrate); - dealii::FEValues fe_values_extra(dg->mapping_collection[poly_degree], dg->fe_collection[poly_degree], quad_extra, + dealii::FEValues fe_values_extra(*(dg->high_order_grid.mapping_fe_field), dg->fe_collection[poly_degree], quad_extra, dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points); const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points; std::array soln_at_q; diff --git a/src/testing/euler_gaussian_bump.cpp b/src/testing/euler_gaussian_bump.cpp index 1bb4ecb20..bd36c9c68 100644 --- a/src/testing/euler_gaussian_bump.cpp +++ b/src/testing/euler_gaussian_bump.cpp @@ -191,7 +191,7 @@ ::run_test () const //n_subdivisions[1] = n_1d_cells[0]; // y-direction //n_subdivisions[0] = 4*n_subdivisions[1]; // x-direction n_subdivisions[1] = n_1d_cells[0]; // y-direction - n_subdivisions[0] = 3*n_subdivisions[1]; // x-direction + n_subdivisions[0] = 9*n_subdivisions[1]; // x-direction dealii::Point<2> p1(-1.5,0.0), p2(1.5,y_height); const bool colorize = true; dealii::parallel::distributed::Triangulation grid(this->mpi_communicator, @@ -224,7 +224,6 @@ ::run_test () const // Create DG object std::shared_ptr < DGBase > dg = DGFactory::create_discontinuous_galerkin(¶m, poly_degree, &grid); - // Initialize coarse grid solution with free-stream dg->allocate_system (); dealii::VectorTools::interpolate(dg->dof_handler, initial_conditions, dg->solution); @@ -240,7 +239,9 @@ ::run_test () const dealii::LinearAlgebra::distributed::Vector old_solution(dg->solution); dealii::parallel::distributed::SolutionTransfer, dealii::hp::DoFHandler> solution_transfer(dg->dof_handler); solution_transfer.prepare_for_coarsening_and_refinement(old_solution); + dg->high_order_grid.prepare_for_coarsening_and_refinement(); grid.refine_global (1); + dg->high_order_grid.execute_coarsening_and_refinement(); dg->allocate_system (); dg->solution.zero_out_ghosts(); solution_transfer.interpolate(dg->solution); @@ -346,12 +347,13 @@ ::run_test () const const double last_slope = log(entropy_error[n_grids-1]/entropy_error[n_grids-2]) / log(grid_size[n_grids-1]/grid_size[n_grids-2]); - double before_last_slope = last_slope; - if ( n_grids > 2 ) { - before_last_slope = log(entropy_error[n_grids-2]/entropy_error[n_grids-3]) - / log(grid_size[n_grids-2]/grid_size[n_grids-3]); - } - const double slope_avg = 0.5*(before_last_slope+last_slope); + //double before_last_slope = last_slope; + //if ( n_grids > 2 ) { + // before_last_slope = log(entropy_error[n_grids-2]/entropy_error[n_grids-3]) + // / log(grid_size[n_grids-2]/grid_size[n_grids-3]); + //} + //const double slope_avg = 0.5*(before_last_slope+last_slope); + const double slope_avg = last_slope; const double slope_diff = slope_avg-expected_slope; double slope_deficit_tolerance = -std::abs(manu_grid_conv_param.slope_deficit_tolerance); diff --git a/src/testing/grid_study.cpp b/src/testing/grid_study.cpp index ecd46844c..66dc2da32 100644 --- a/src/testing/grid_study.cpp +++ b/src/testing/grid_study.cpp @@ -63,7 +63,7 @@ ::integrate_solution_over_domain(DGBase &dg) const int overintegrate = 10; dealii::QGauss quad_extra(dg.max_degree+1+overintegrate); //dealii::MappingQ mappingq_temp(dg.max_degree+1); - dealii::FEValues fe_values_extra(dg.mapping_collection[dg.max_degree], dg.fe_collection[dg.max_degree], quad_extra, + dealii::FEValues fe_values_extra(*(dg.high_order_grid.mapping_fe_field), dg.fe_collection[dg.max_degree], quad_extra, dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points); const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points; std::array soln_at_q; @@ -260,7 +260,7 @@ ::run_test () const int overintegrate = 10; dealii::QGauss quad_extra(dg->max_degree+1+overintegrate); //dealii::MappingQ mappingq(dg->max_degree+1); - dealii::FEValues fe_values_extra(dg->mapping_collection[poly_degree], dg->fe_collection[poly_degree], quad_extra, + dealii::FEValues fe_values_extra(*(dg->high_order_grid.mapping_fe_field), dg->fe_collection[poly_degree], quad_extra, dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points); const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points; std::array soln_at_q; diff --git a/tests/integration_tests_control_files/euler_integration/2d_euler_gaussian_bump.prm b/tests/integration_tests_control_files/euler_integration/2d_euler_gaussian_bump.prm index ca7c8a50f..6b8970ea5 100644 --- a/tests/integration_tests_control_files/euler_integration/2d_euler_gaussian_bump.prm +++ b/tests/integration_tests_control_files/euler_integration/2d_euler_gaussian_bump.prm @@ -60,9 +60,9 @@ subsection manufactured solution convergence study set grid_progression_add = 0 # Initial grid of size (initial_grid_size)^dim - set initial_grid_size = 3 + set initial_grid_size = 7 # Number of grids in grid study - set number_of_grids = 4 + set number_of_grids = 3 end diff --git a/tests/unit_tests/grid/high_order_grid_test.cpp b/tests/unit_tests/grid/high_order_grid_test.cpp index 5d7000335..53323ff8a 100644 --- a/tests/unit_tests/grid/high_order_grid_test.cpp +++ b/tests/unit_tests/grid/high_order_grid_test.cpp @@ -94,7 +94,8 @@ int main (int argc, char * argv[]) HighOrderGrid high_order_grid(&all_parameters, poly_degree, &grid); - dealii::MappingFEField, dealii::DoFHandler> mapping = high_order_grid.get_MappingFEField(); + //std::shared_ptr < dealii::MappingFEField, dealii::DoFHandler> > mapping = (high_order_grid.mapping_fe_field); + auto mapping = (high_order_grid.mapping_fe_field); grid.reset_all_manifolds(); dealii::ConvergenceTable convergence_table; @@ -138,7 +139,7 @@ int main (int argc, char * argv[]) data_out.add_data_vector(high_order_grid.nodes, solution_names); const int iproc = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); - data_out.build_patches(mapping, poly_degree, dealii::DataOut::CurvedCellRegion::curved_inner_cells); + data_out.build_patches(*mapping, poly_degree, dealii::DataOut::CurvedCellRegion::curved_inner_cells); std::string filename = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D-"; filename += "Degree" + dealii::Utilities::int_to_string(poly_degree, 2) + "."; filename += dealii::Utilities::int_to_string(igrid, 4) + "."; @@ -197,7 +198,7 @@ int main (int argc, char * argv[]) const unsigned int n_quad_pts_1D = poly_degree+1+overintegrate; dealii::QGauss quadrature(n_quad_pts_1D); const unsigned int n_quad_pts = quadrature.size(); - dealii::FEValues fe_values(mapping, high_order_grid.fe_system, quadrature, + dealii::FEValues fe_values(*mapping, high_order_grid.fe_system, quadrature, dealii::update_jacobians | dealii::update_JxW_values | dealii::update_quadrature_points);