From 767fe8b34ebc756eaea6e67e479aa4ba05757c35 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 14 Mar 2018 00:45:22 +0100 Subject: [PATCH] Fix CeresFE --- CeresFE/src/ceres.cc | 33 ++++++++++++++-------------- CeresFE/support_code/ellipsoid_fit.h | 4 +--- CeresFE/support_code/local_math.h | 2 -- 3 files changed, 17 insertions(+), 22 deletions(-) diff --git a/CeresFE/src/ceres.cc b/CeresFE/src/ceres.cc index d47a70e..914cc4e 100644 --- a/CeresFE/src/ceres.cc +++ b/CeresFE/src/ceres.cc @@ -185,7 +185,7 @@ private: BlockVector solution; BlockVector system_rhs; - std_cxx1x::shared_ptr::type> A_preconditioner; + std::shared_ptr::type> A_preconditioner; ellipsoid_fit ellipsoid; }; @@ -228,6 +228,10 @@ class RightHandSide: public Function { public: RightHandSide () : Function(dim+1) {} + using Function::value; + using Function::vector_value; + using Function::vector_value_list; + virtual double value(const Point &p, const unsigned int component, A_Grav_namespace::AnalyticGravity *aGrav) const; @@ -352,8 +356,8 @@ void SchurComplement::vmult(Vector &dst, template StokesProblem::StokesProblem(const unsigned int degree) : degree(degree), - mapping(), triangulation(Triangulation::maximum_smoothing), + mapping(), fe(FE_Q(degree + 1), dim, FE_Q(degree), 1), dof_handler(triangulation), quadrature_formula(degree + 2), @@ -384,7 +388,7 @@ void StokesProblem::setup_dofs() { BoundaryValuesP(), constraints, component_maskP); } { - std::set no_normal_flux_boundaries; + std::set no_normal_flux_boundaries; no_normal_flux_boundaries.insert(99); VectorTools::compute_no_normal_flux_constraints(dof_handler, 0, no_normal_flux_boundaries, constraints); @@ -392,7 +396,7 @@ void StokesProblem::setup_dofs() { constraints.close(); - std::vector dofs_per_block(2); + std::vector dofs_per_block(2); DoFTools::count_dofs_per_block(dof_handler, dofs_per_block, block_component); n_u = dofs_per_block[0]; @@ -404,7 +408,7 @@ void StokesProblem::setup_dofs() { << std::endl; { - BlockCompressedSimpleSparsityPattern csp(2, 2); + BlockDynamicSparsityPattern csp(2, 2); csp.block(0, 0).reinit(n_u, n_u); csp.block(1, 0).reinit(n_p, n_u); @@ -464,7 +468,6 @@ double Rheology::get_eta(double &r, double &z) double ecc = system_parameters::q_axes[0] / system_parameters::p_axes[0]; double Rminusr = system_parameters::q_axes[0] - system_parameters::p_axes[0]; double approx_a = std::sqrt(r * r + z * z * ecc * ecc); - double approx_b = approx_a / ecc; double group1 = r * r + z * z - Rminusr * Rminusr; double a0 = approx_a; @@ -600,7 +603,6 @@ void StokesProblem::initialize_eta_and_G() { //defines local viscosity double local_viscosity = 0; - unsigned int m_id = cell->material_id(); local_viscosity = rheology.get_eta(r_value, z_value); @@ -639,7 +641,7 @@ void StokesProblem::assemble_system() { FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); Vector local_rhs(dofs_per_cell); - std::vector local_dof_indices(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); // runs the gravity script function const RightHandSide right_hand_side; @@ -741,7 +743,6 @@ void StokesProblem::assemble_system() { double &local_old_phiphi_stress = local_quadrature_points_history[q].old_phiphi_stress; double r_value = fe_values.quadrature_point(q)[0]; - double z_value = fe_values.quadrature_point(q)[1]; // get local density based on mat id double local_density = system_parameters::rho[m_id]; @@ -836,7 +837,6 @@ void StokesProblem::assemble_system() { double &local_old_phiphi_stress = local_quadrature_points_history[q].old_phiphi_stress; double r_value = fe_values.quadrature_point(q)[0]; - double z_value = fe_values.quadrature_point(q)[1]; double local_density = system_parameters::rho[m_id]; @@ -930,7 +930,7 @@ void StokesProblem::assemble_system() { std::cout << " Computing preconditioner..." << std::endl << std::flush; - A_preconditioner = std_cxx1x::shared_ptr< + A_preconditioner = std::shared_ptr< typename InnerPreconditioner::type>( new typename InnerPreconditioner::type()); A_preconditioner->initialize(system_matrix.block(0, 0), @@ -1239,7 +1239,6 @@ void StokesProblem::solution_stesses() { } // Finds adjusted effective viscosity - double cell_effective_viscosity = 0; if (system_parameters::plasticity_on) { if(system_parameters::failure_criterion == 0) //Apply Byerlee's rule @@ -1433,7 +1432,7 @@ void StokesProblem::smooth_eta_field(std::vector failing_cells) std::vector cell_touched(triangulation.n_active_cells(), false); std::vector< TriaIterator< CellAccessor > > neighbor_cells; std::vector neighbor_indices; - int start_cell = 0; // Which cell in the neighbor_cells vector to start from + unsigned int start_cell = 0; // Which cell in the neighbor_cells vector to start from int new_cells_found = 0; neighbor_cells.push_back(cell); neighbor_indices.push_back(cell_no); @@ -1441,7 +1440,7 @@ void StokesProblem::smooth_eta_field(std::vector failing_cells) while(find_more_cells) { new_cells_found = 0; - for(int i = start_cell; i::faces_per_cell; ++f) { @@ -1741,7 +1740,7 @@ void StokesProblem::write_vertices(unsigned char boundary_that_we_need) { triangulation.begin_active(); cell != triangulation.end(); ++cell) for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) { - unsigned char boundary_ids = cell->face(f)->boundary_indicator(); + unsigned char boundary_ids = cell->face(f)->boundary_id(); if(boundary_ids == boundary_that_we_need) { for (unsigned int v=0; v::vertices_per_face; ++v) @@ -1834,9 +1833,9 @@ void StokesProblem::setup_initial_mesh() { if (fabs(cell->face(f)->center()[i]) > zero_tolerance) how_many++; if (how_many==dim) - cell->face(f)->set_all_boundary_indicators(0); // if face center coordinates > zero_tol, set bnry indicators to 0 + cell->face(f)->set_all_boundary_ids(0); // if face center coordinates > zero_tol, set bnry indicators to 0 else - cell->face(f)->set_all_boundary_indicators(99); + cell->face(f)->set_all_boundary_ids(99); } } } diff --git a/CeresFE/support_code/ellipsoid_fit.h b/CeresFE/support_code/ellipsoid_fit.h index d47f3e3..c99458b 100644 --- a/CeresFE/support_code/ellipsoid_fit.h +++ b/CeresFE/support_code/ellipsoid_fit.h @@ -11,7 +11,6 @@ #include #include #include -#include #include #include #include @@ -80,14 +79,13 @@ void ellipsoid_fit::compute_fit(std::vector &ell, unsigned char bou std::vector ind_bnry_col; // assemble the sensitivity matrix and r.h.s. - double zero_tolerance = 1e-3; for (; cell != endc; ++cell) { if (boundary_that_we_need != 0) cell->set_manifold_id(cell->material_id()); for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) { if (boundary_that_we_need == 0) //if this is the outer surface, then look for boundary ID 0; otherwise look for material ID change. { - boundary_ids = cell->face(f)->boundary_indicator(); + boundary_ids = cell->face(f)->boundary_id(); if (boundary_ids == boundary_that_we_need) { for (unsigned int v = 0; v < GeometryInfo::vertices_per_face; ++v) diff --git a/CeresFE/support_code/local_math.h b/CeresFE/support_code/local_math.h index 2dcc5a9..cbece34 100644 --- a/CeresFE/support_code/local_math.h +++ b/CeresFE/support_code/local_math.h @@ -11,8 +11,6 @@ #define PI 3.14159265358979323846 #define TWOPI 6.283185307179586476925287 #define SECSINYEAR 3.155692608e+07 -#define MAX(x,y) ((x) > (y)) ? (x) : (y) -#define MIN(x,y) ((x) < (y)) ? (x) : (y) //#define ABS(a) ((a) < 0 ? -(a) : (a)) //double factorial(int n) -- 2.39.5