{"action":"closed","issue":{"url":"https://api.github.com/repos/dealii/dealii/issues/2469","repository_url":"https://api.github.com/repos/dealii/dealii","labels_url":"https://api.github.com/repos/dealii/dealii/issues/2469/labels{/name}","comments_url":"https://api.github.com/repos/dealii/dealii/issues/2469/comments","events_url":"https://api.github.com/repos/dealii/dealii/issues/2469/events","html_url":"https://github.com/dealii/dealii/issues/2469","id":146311226,"number":2469,"title":"Error when plotting curved interior cells with vtk","user":{"login":"bangerth","id":7504421,"avatar_url":"https://avatars.githubusercontent.com/u/7504421?v=3","gravatar_id":"","url":"https://api.github.com/users/bangerth","html_url":"https://github.com/bangerth","followers_url":"https://api.github.com/users/bangerth/followers","following_url":"https://api.github.com/users/bangerth/following{/other_user}","gists_url":"https://api.github.com/users/bangerth/gists{/gist_id}","starred_url":"https://api.github.com/users/bangerth/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bangerth/subscriptions","organizations_url":"https://api.github.com/users/bangerth/orgs","repos_url":"https://api.github.com/users/bangerth/repos","events_url":"https://api.github.com/users/bangerth/events{/privacy}","received_events_url":"https://api.github.com/users/bangerth/received_events","type":"User","site_admin":false},"labels":[{"id":114738732,"url":"https://api.github.com/repos/dealii/dealii/labels/Bug","name":"Bug","color":"fc2929","default":false}],"state":"closed","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":5,"created_at":"2016-04-06T13:17:59Z","updated_at":"2016-11-16T00:54:22Z","closed_at":"2016-11-16T00:54:22Z","body":"As reported here, https://groups.google.com/forum/#!topic/dealii/_uk1JDakBLc , there is a bug somewhere in `DataOut::build_patches` or `DataOut::build_one_patch` when using curved geometries. My best guess is that we compute the location of points on the edges of the inner square using the manifold of the square, rather than the manifold of the edge. This then leads us to assume straight edges for the inner square but curved edges for the outer cells.\r\n\r\nTest program by Sebastian Gonzalez-Pintor:\r\n\r\n```\r\n\r\n// Testing grids with Manifolds by approximating pi.\r\n\r\n#include <deal.II/base/quadrature_lib.h>\r\n#include <deal.II/base/convergence_table.h>\r\n\r\n#include <deal.II/grid/grid_generator.h>\r\n#include <deal.II/grid/tria_accessor.h>\r\n#include <deal.II/grid/tria_iterator.h>\r\n#include <deal.II/grid/tria_boundary_lib.h>\r\n#include <deal.II/grid/tria.h>\r\n#include <deal.II/grid/manifold_lib.h>\r\n#include <deal.II/grid/grid_out.h>\r\n\r\n#include <deal.II/dofs/dof_handler.h>\r\n#include <deal.II/dofs/dof_accessor.h>\r\n\r\n#include <deal.II/fe/fe_q.h>\r\n#include <deal.II/fe/fe_values.h>\r\n#include <deal.II/fe/mapping_q.h>\r\n\r\n#include <deal.II/numerics/data_out.h>\r\n\r\n// And this again is C++:\r\n#include <iostream>\r\n#include <fstream>\r\n#include <cmath>\r\n\r\nusing namespace dealii;\r\nusing namespace std;\r\n\r\n// mainly Step-10 without computing the perimeter\r\nnamespace Step10\r\n{\r\n  using namespace dealii;\r\n\r\n  const long double pi = 3.141592653589793238462643;\r\n  const double inner_radius = 1./(1. + std::sqrt(2));\r\n  const double true_inner_area = pi*std::pow(inner_radius,2);\r\n  const double true_outer_area = pi - true_inner_area;\r\n\r\n  void compute_pi_by_area (const unsigned int n_degree)\r\n  {\r\n    std::cout << \"Computation of Pi by the area:\" << std::endl\r\n              << \"==============================\" << std::endl;\r\n    printf(\"[deg] \\t[ TArea ] \\t[TA - PI] \\t[err(OA)] \\t[err(IA)]\\n\\n\");\r\n\r\n    Vector<double> zeros;\r\n\r\n    const Point<2> center (0,0);\r\n    const SphericalManifold<2> boundary_description (center);\r\n\r\n    for (unsigned int degree=1; degree <= n_degree; ++degree)\r\n    {\r\n        const QGauss<2> quadrature(degree);\r\n        Triangulation<2> triangulation;\r\n        GridGenerator::hyper_ball (triangulation, Point< 2 >(0.0,0.0), 1.0 );\r\n\r\n        const MappingQ<2> mapping (degree, true);\r\n        const FE_Q<2>     dummy_fe (1);\r\n        DoFHandler<2> dof_handler (triangulation);\r\n        FEValues<2> fe_values (mapping, dummy_fe, quadrature,\r\n                             update_quadrature_points | update_JxW_values);\r\n\r\n        for (Triangulation<2>::active_cell_iterator\r\n          cell = triangulation.begin_active(); cell!=triangulation.end(); ++cell) \r\n          {\r\n            if (cell->center().square() < 1.e-5)\r\n              cell->set_material_id(1);\r\n            else\r\n              cell->set_material_id(0);\r\n          }\r\n\r\n        Triangulation<2>::active_cell_iterator\r\n          cell = triangulation.begin_active(),\r\n          endc = triangulation.end();\r\n\r\n        for (; cell!=endc; ++cell) \r\n          if (cell->material_id() != 1)\r\n            cell->set_all_manifold_ids (100);\r\n\r\n        triangulation.set_manifold (100, boundary_description);\r\n\r\n        dof_handler.distribute_dofs (dummy_fe); \r\n\r\n        zeros.reinit(triangulation.n_active_cells());\r\n        cell = dof_handler.begin_active();\r\n        for (unsigned int i = 0; cell!=endc; ++cell, ++i)\r\n        {\r\n          zeros[i] = static_cast<double>(cell->material_id());\r\n        }\r\n\r\n        double area_inner = 0;\r\n        double area_outer = 0;\r\n        cell = dof_handler.begin_active();\r\n        for (; cell!=endc; ++cell)\r\n        {\r\n          fe_values.reinit (cell);\r\n          for (unsigned int i=0; i < fe_values.n_quadrature_points; ++i)\r\n          {\r\n            if (cell->material_id() == 1)\r\n              area_inner += fe_values.JxW (i);\r\n            else\r\n            area_outer += fe_values.JxW (i);\r\n          }\r\n        }\r\n        const double area_total = area_outer + area_inner;\r\n\r\n        printf(\" %d \\t%.3e \\t%.3e \\t%.3e \\t%.3e\\n\",\r\n               degree,\r\n               (double)area_total,\r\n               (double)std::fabs(pi - area_total),\r\n               (double)(area_outer - true_outer_area),\r\n               (double)(area_inner - true_inner_area));\r\n\r\n        { // plotting the grid\r\n          DataOut<2> data_out;\r\n          data_out.attach_dof_handler (dof_handler);\r\n          std::ostringstream ss;              \r\n          ss <<  \"d2\" << \"_p\" << degree;\r\n          data_out.add_data_vector (zeros, std::string(\"grid\") );\r\n          std::ofstream o_vtk ( ss.str() + \".vtk\" );\r\n          data_out.build_patches (mapping, 5, DataOut<2>::curved_inner_cells );\r\n          data_out.write_vtk (o_vtk);\r\n        }\r\n    } // end for degree\r\n  printf(\"\\n\");\r\n  } // end compute_pi_by_area\r\n\r\n} // end namespace\r\n\r\nint main ()\r\n{\r\n  try\r\n    {\r\n      std::cout.precision (16);\r\n\r\n      const unsigned int n_degree = 4;\r\n\r\n      Step10::compute_pi_by_area (n_degree);\r\n    }\r\n  catch (std::exception &exc)\r\n    {\r\n      std::cerr << std::endl << std::endl\r\n                << \"----------------------------------------------------\"\r\n                << std::endl;\r\n      std::cerr << \"Exception on processing: \" << std::endl\r\n                << exc.what() << std::endl\r\n                << \"Aborting!\" << std::endl\r\n                << \"----------------------------------------------------\"\r\n                << std::endl;\r\n\r\n      return 1;\r\n    }\r\n  catch (...)\r\n    {\r\n      std::cerr << std::endl << std::endl\r\n                << \"----------------------------------------------------\"\r\n                << std::endl;\r\n      std::cerr << \"Unknown exception!\" << std::endl\r\n                << \"Aborting!\" << std::endl\r\n                << \"----------------------------------------------------\"\r\n                << std::endl;\r\n      return 1;\r\n    }\r\n  return 0;\r\n}\r\n```\r\n"}}
{"actor":{"display_login":"masterleinad"}}