From cfed2b3a858307755dcd339a649e898ec4beea98 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 23 Apr 2002 08:44:12 +0000 Subject: [PATCH] More test. Fixes. git-svn-id: https://svn.dealii.org/trunk@5713 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-14/step-14.cc | 208 ++++++++++++++-------------- 1 file changed, 101 insertions(+), 107 deletions(-) diff --git a/deal.II/examples/step-14/step-14.cc b/deal.II/examples/step-14/step-14.cc index 930a01253f..a821c34623 100644 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@ -1,3 +1,5 @@ +// TODO: bv for primal problem wrong! +// check Galerkin orthogonality /* $Id$ */ /* Author: Wolfgang Bangerth, ETH Zurich, 2002 */ @@ -149,8 +151,8 @@ namespace Evaluation results_table.add_value ("u(x_0)", point_value); std::cout << " Point value=" << point_value - << ", exact value=1, error=" - << 1.-point_value << std::endl; + << ", exact value=1.59492, error=" + << 1.594915543-point_value << std::endl; }; }; @@ -449,7 +451,7 @@ namespace LaplaceSolver void Solver::LinearSystem::solve (Vector &solution) const { - SolverControl solver_control (1000, 1e-12); + SolverControl solver_control (1000, 1e-16); PrimitiveVectorMemory<> vector_memory; SolverCG<> cg (solver_control, vector_memory); @@ -723,28 +725,28 @@ double RightHandSide::value (const Point &p, const unsigned int /*component*/) const { -// double q = p(0); -// for (unsigned int i=1; i::active_cell_iterator> FaceEntry; - typedef typename std::map::face_iterator,FaceEntry> FaceIntegrals; + typedef + typename std::map::face_iterator,double> + FaceIntegrals; /** @@ -1167,13 +1180,13 @@ namespace LaplaceSolver const Quadrature &quadrature, const Quadrature &face_quadrature, const Function &rhs_function, - const Function &boundary_values, + const Function &bv, const DualFunctional::DualFunctionalBase &dual_functional) : Base (coarse_grid), PrimalSolver (coarse_grid, primal_fe, quadrature, face_quadrature, - rhs_function, boundary_values), + rhs_function, bv), DualSolver (coarse_grid, dual_fe, quadrature, face_quadrature, dual_functional) @@ -1212,10 +1225,6 @@ namespace LaplaceSolver { Vector error_indicators (triangulation->n_active_cells()); estimate_error (error_indicators); - std::cout << " Estimated error=" - << std::accumulate (error_indicators.begin(), - error_indicators.end(), 0.) - << std::endl; DataOut data_out; ofstream x("x"); Vector xe (error_indicators.begin(), @@ -1223,7 +1232,7 @@ namespace LaplaceSolver data_out.attach_dof_handler (DualSolver::dof_handler); data_out.add_data_vector (xe, "e"); data_out.build_patches (); - data_out.write_gnuplot (x); + data_out.write_gmv (x); std::transform (error_indicators.begin(), error_indicators.end(), @@ -1236,50 +1245,42 @@ namespace LaplaceSolver }; + template void WeightedResidual::output_solution () const { - for (unsigned int solution_type=0; solution_type<2; ++solution_type) - { - DataOut data_out; + Vector primal_solution (DualSolver::dof_handler.n_dofs()); + FETools::interpolate (PrimalSolver::dof_handler, + PrimalSolver::solution, + DualSolver::dof_handler, + primal_solution); - switch (solution_type) - { - case 0: - data_out.attach_dof_handler (PrimalSolver::dof_handler); - data_out.add_data_vector (PrimalSolver::solution, - "primal_solution"); - break; - case 1: - data_out.attach_dof_handler (DualSolver::dof_handler); - data_out.add_data_vector (DualSolver::solution, - "dual_solution"); - break; - default: - Assert (false, ExcInternalError()); - }; - data_out.build_patches (); + DataOut data_out; + data_out.attach_dof_handler (DualSolver::dof_handler); + data_out.add_data_vector (primal_solution, + "primal_solution"); + data_out.add_data_vector (DualSolver::solution, + "dual_solution"); + + data_out.build_patches (1); #ifdef HAVE_STD_STRINGSTREAM - std::ostringstream filename; + std::ostringstream filename; #else - std::ostrstream filename; + std::ostrstream filename; #endif - filename << "solution-" - << (solution_type == 0 ? - "primal-" : "dual-") - << refinement_cycle - << ".gnuplot" - << std::ends; + filename << "solution-" + << refinement_cycle + << ".gnuplot" + << std::ends; #ifdef HAVE_STD_STRINGSTREAM - std::ofstream out (filename.str().c_str()); + std::ofstream out (filename.str().c_str()); #else - std::ofstream out (filename.str()); + std::ofstream out (filename.str()); #endif - data_out.write (out, DataOut::gnuplot); - }; + data_out.write (out, DataOut::gnuplot); }; @@ -1365,7 +1366,7 @@ namespace LaplaceSolver for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - face_integrals[cell->face(face_no)].first = -1e20; + face_integrals[cell->face(face_no)] = -1e20; // Then set up a vector with // error indicators. Reserve one @@ -1399,7 +1400,7 @@ namespace LaplaceSolver // this, note that the cell terms // are already set, and that only // the edge terms need to be - // collected. For this, loop over + // collected. Thus, loop over // all cells and their faces, // make sure that the // contributions of each of the @@ -1412,17 +1413,16 @@ namespace LaplaceSolver for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { - Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(), + Assert(face_integrals.find(cell->face(face_no)) != + face_integrals.end(), ExcInternalError()); - if (true || (face_integrals[cell->face(face_no)].second - == - cell)) - error_indicators(present_cell) - += -0.5*face_integrals[cell->face(face_no)].first; - else - error_indicators(present_cell) - -= -0.5*face_integrals[cell->face(face_no)].first; + error_indicators(present_cell) + += 0.5*face_integrals[cell->face(face_no)]; }; + std::cout << " Estimated error=" + << std::accumulate (error_indicators.begin(), + error_indicators.end(), 0.) + << std::endl; }; @@ -1536,24 +1536,20 @@ namespace LaplaceSolver // enter this face into // the list of faces with // a zero contribution to - // the error, and also - // mark the cell on which - // we computed this - // value. + // the error. if (cell->face(face_no)->at_boundary()) { - face_integrals[cell->face(face_no)].first = 0; - face_integrals[cell->face(face_no)].second = cell; + face_integrals[cell->face(face_no)] = 0; continue; }; // Next, note that since // we want to compute the // jump terms on each - // face only once, + // face only once // although we access it - // twice if it is not at - // the boundary, we have + // twice (if it is not at + // the boundary), we have // to define some rules // who is responsible for // computing on a face: @@ -1635,10 +1631,10 @@ namespace LaplaceSolver // next cell for this // thread. Note again that // the cells for each of the - // threads are - // interleaved. If we are at - // the end of our workload, - // jump out of the loop. + // threads are interleaved. + // If we are at the end of + // our workload, jump out + // of the loop. for (unsigned int t=0; ((t::dof_handler.end())); ++t, ++cell, ++cell_index); @@ -1682,7 +1678,9 @@ namespace LaplaceSolver cell_data.dual_weights); // ...and finally build the sum - // over all quadrature points: + // over all quadrature points and + // store it with the present + // cell: double sum = 0; for (unsigned int p=0; pface(face_no)) != face_integrals.end(), ExcInternalError()); - Assert (face_integrals[cell->face(face_no)].first == -1e20, + Assert (face_integrals[cell->face(face_no)] == -1e20, ExcInternalError()); // ...then store computed value // at assigned location: - face_integrals[cell->face(face_no)].first = face_integral; - face_integrals[cell->face(face_no)].second = cell; + face_integrals[cell->face(face_no)] = face_integral; }; @@ -1938,10 +1935,8 @@ namespace LaplaceSolver face_integral += (face_data.jump_residual[p] * face_data.dual_weights[p] * face_data.fe_face_values_neighbor.JxW(p)); - face_integrals[neighbor_child->face(neighbor_neighbor)].first + face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral; - face_integrals[neighbor_child->face(neighbor_neighbor)].second - = cell; }; // Once the contributions of all @@ -1964,15 +1959,14 @@ namespace LaplaceSolver Assert (face_integrals.find(face->child(subface_no)) != face_integrals.end(), ExcInternalError()); - Assert (face_integrals[face->child(subface_no)].first != -1e20, + Assert (face_integrals[face->child(subface_no)] != -1e20, ExcInternalError()); - sum += face_integrals[face->child(subface_no)].first; + sum += face_integrals[face->child(subface_no)]; }; // Finally store the value with // the parent face. - face_integrals[face].first = sum; - face_integrals[face].second = cell; + face_integrals[face] = sum; }; }; @@ -2006,7 +2000,7 @@ run_simulation (LaplaceSolver::Base &solver, }; - if (solver.n_dofs() < 2000) + if (solver.n_dofs() < 5000) solver.refine_grid (); else break; @@ -2074,7 +2068,7 @@ void solve_problem (const std::string &solver_name) Triangulation triangulation; // create_triangulation (triangulation); GridGenerator::hyper_cube (triangulation, -1, 1); - triangulation.refine_global (2); + triangulation.refine_global (5); const FE_Q primal_fe(1); const FE_Q dual_fe(2); const QGauss4 quadrature; -- 2.39.5