<https://lh3.googleusercontent.com/-jdyS6vcktFY/V6kDXk0vxGI/AAAAAAAAFR8/CnkWKXhHuE4IMYxgyhO9RxCQ4owNt-F-wCLcB/s1600/ellipsoidal_mesh0001.png> I am trying to obtain the mesh for a codimension-1 ellipsoid and attach an ellipsoidal manifold to it in order to refine the mesh. My strategy is failing and I have a few questions.
*My strategy: * some definitions: dim=2; spacedim=3; chartdim=spacedim-1; a=1; b=2;c=3; Set up codimension-1 Triangulation: Triangulation<dim,spacedim> tria; Use GridGenerator::hyper_sphere() to obtain a codimension 1 mesh: GridGenerator::hyper_sphere (tria, Point<spacedim>(0.0,0.0,0.0), 1.0); Use GridTools::transform() to map triangulation of sphere onto triangulation of ellipsoid (I have a question about this, see below). Here, grid_transform has signature Point<spacedim> grid_transform(const Point<spacedim> X) GridTools::transform(&grid_transform, tria); Next, subclass ChartManifold<dim,spacedim,chartdim> in order to obtain a push_forward() and a pull_back() function for my parametric ellipsoidal manifold description. Then attach that manifold to the triangulation: Ellipsoid<dim,spacedim,chartdim> ellipsoid(a,b,c); tria.set_manifold(0,ellipsoid); Now refine the triangulation in order to see a refined ellipsoidal mesh: tria.refine_global(1); *The Result: *Complete garbage! I also don't really have an idea why this is not working - I've checked the dimensions on my push_forward() and pull_back() functions, I think that I'm setting the periodicity correctly in my ChartManifold superclass, and I've also checked that my simple math is okay. Let me know what other information would be useful. Below is the refined mesh for a=b=c=1. This should just be a sphere. <https://lh3.googleusercontent.com/-jdyS6vcktFY/V6kDXk0vxGI/AAAAAAAAFR8/CnkWKXhHuE4IMYxgyhO9RxCQ4owNt-F-wCLcB/s1600/ellipsoidal_mesh0001.png> *The code:* template <int dim,int spacedim,int chartdim=spacedim-1> class Ellipsoid: public ChartManifold<dim,spacedim,chartdim> { public: Ellipsoid(double,double,double); Point<chartdim> pull_back(const Point<spacedim> &space_point) const; Point<spacedim> push_forward(const Point<chartdim> &chart_point) const; private: double a,b,c; double max_axis; const Point<spacedim> center; SphericalManifold<dim,spacedim> reference_sphere; }; template <int dim, int spacedim, int chartdim> Ellipsoid<dim,spacedim,chartdim>::Ellipsoid(double a, double b, double c): ChartManifold<dim,spacedim,chartdim>(Point<chartdim>(2*numbers::PI, 2*numbers::PI)), a(a), b(b),c(c), center(0,0,0), reference_sphere(center) { max_axis = std::max(std::max(a,b),c); } template <int dim,int spacedim, int chartdim> Point<chartdim> Ellipsoid<dim,spacedim,chartdim>::pull_back(const Point<spacedim> &space_point) const { double x,y,z, u,v,w; // get point on ellipsoid x = space_point[0]; y = space_point[1]; z = space_point[2]; std::cout << "using a,b,c: " << std::endl; std::cout << a << " " << b << " " << c << std::endl; std::cout << "from pull_back: " << std::endl; std::cout << "space_point: " << std::endl; std::cout << x << " " << y << " " << z << std::endl; // map ellipsoid point onto sphere u = x/a; v = y/b; w = z/c; std::cout << "pulls back to : " << std::endl; std::cout << u << " " << v << " " << w << std::endl; std::cout << "on sphere." << std::endl; Point<spacedim> p(u,v,w); // use reference_sphere's pull_back function Point<spacedim> q = reference_sphere.pull_back(p); Point<chartdim> chart_point; std::cout << "sphere pull_back: " << std::endl; std::cout << q[0] << " " << q[1] << " " << q[2] << std::endl; std::cout << "r theta phi" << std::endl; std::cout << "..........." << std::endl; chart_point[0] = q[1]; chart_point[1] = q[2]; // return (theta,phi) in the chart domain return chart_point; } template <int dim,int spacedim, int chartdim> Point<spacedim> Ellipsoid<dim,spacedim,chartdim>::push_forward(const Point<chartdim> &chart_point) const { double theta,phi, x,y,z; phi = chart_point[0]; theta = chart_point[1]; Point<spacedim> p(max_axis,theta,phi); // map theta,phi in chart domain onto reference_sphere with radius max_axis Point<spacedim> X = reference_sphere.push_forward(p); // map point on sphere onto ellipsoid x = a*X[0]; y = b*X[1]; z = c*X[2]; Point<spacedim> space_point(x,y,z); // return point on ellipsoid return space_point; } Point<3> grid_transform (const Point<3> &X) { // transform points on sphere onto ellipsoid double a,b,c; a = 1.0; b = 1.0; c = 1.0; double x,y,z; x = a*X(0); y = b*X(1); z = c*X(2); return Point<3>(x,y,z); } void assemble_mesh_and_manifold() { const int dim = 2; const int spacedim = 3; const int chartdim = 2; double a,b,c; a = 1; b=1; c=1; Ellipsoid<dim,spacedim,chartdim> ellipsoid(a,b,c); Triangulation<dim,spacedim> tria; // generate coarse spherical mesh GridGenerator::hyper_sphere (tria, Point<spacedim>(0.0,0.0,0.0), 1.0); for (Triangulation<dim,spacedim>::active_cell_iterator cell=tria.begin_active(); cell!=tria.end(); ++cell) cell->set_all_manifold_ids(0); print_mesh_info(tria, "spherical_mesh.vtk"); GridTools::transform(&grid_transform, tria); tria.set_manifold(0,ellipsoid); tria.refine_global(1); print_mesh_info(tria, "ellipsoidal_mesh.vtk"); } int main () { assemble_mesh_and_manifold(); } Attached is the .vtk of my refined ellipsoidal mesh and the .cc file that generates this output (mostly reproduced above. Thank you, Tom -- 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.
/* ---------------------------------------------------------------------
*
* sphere_to_ellipsoid.cc_
modified from step-49 of the tutorial July 25, 2016
Purpose: to generate a grid using dealii tools, modify it according to a
smooth function, and then attach a manifold to it.
* Copyright (C) 2013 - 2015 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Timo Heister, Texas A&M University, 2013
*/
// This tutorial program is odd in the sense that, unlike for most other
// steps, the introduction already provides most of the information on how to
// use the various strategies to generate meshes. Consequently, there is
// little that remains to be commented on here, and we intersperse the code
// with relatively little text. In essence, the code here simply provides a
// reference implementation of what has already been described in the
// introduction.
// @sect3{Include files}
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <iostream>
#include <fstream>
#include <map>
using namespace dealii;
// @sect3{Generating output for a given mesh}
// The following function generates some output for any of the meshes we will
// be generating in the remainder of this program. In particular, it generates
// the following information:
//
// - Some general information about the number of space dimensions in which
// this mesh lives and its number of cells.
// - The number of boundary faces that use each boundary indicator, so that
// it can be compared with what we expect.
//
// Finally, the function outputs the mesh in encapsulated postscript (EPS)
// format that can easily be visualized in the same way as was done in step-1.
template <int dim,int spacedim>
void print_mesh_info(const Triangulation<dim,spacedim> &tria,
const std::string &filename)
{
std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << tria.n_active_cells() << std::endl;
// Next loop over all faces of all cells and find how often each
// boundary indicator is used (recall that if you access an element
// of a std::map object that doesn't exist, it is implicitly created
// and default initialized -- to zero, in the current case -- before
// we then increment it):
{
std::map<unsigned int, unsigned int> boundary_count;
typename Triangulation<dim,spacedim>::active_cell_iterator
cell = tria.begin_active(),
endc = tria.end();
for (; cell!=endc; ++cell)
{
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
{
if (cell->face(face)->at_boundary())
boundary_count[cell->face(face)->boundary_id()]++;
}
}
std::cout << " boundary indicators: ";
for (std::map<unsigned int, unsigned int>::iterator it=boundary_count.begin();
it!=boundary_count.end();
++it)
{
std::cout << it->first << "(" << it->second << " times) ";
}
std::cout << std::endl;
}
// Finally, produce a graphical representation of the mesh to an output
// file:
std::ofstream out (filename.c_str());
GridOut grid_out;
grid_out.write_vtk (tria, out);
std::cout << " written to " << filename
<< std::endl
<< std::endl;
}
template <int dim,int spacedim,int chartdim=spacedim-1>
class Ellipsoid: public ChartManifold<dim,spacedim,chartdim>
{
public:
Ellipsoid(double,double,double);
Point<chartdim> pull_back(const Point<spacedim> &space_point) const;
Point<spacedim> push_forward(const Point<chartdim> &chart_point) const;
private:
double a,b,c;
double max_axis;
const Point<spacedim> center;
SphericalManifold<dim,spacedim> reference_sphere;
};
template <int dim, int spacedim, int chartdim>
Ellipsoid<dim,spacedim,chartdim>::Ellipsoid(double a, double b, double c): ChartManifold<dim,spacedim,chartdim>(Point<chartdim>(2*numbers::PI, 2*numbers::PI)), a(a), b(b),c(c), center(0,0,0), reference_sphere(center)
{
max_axis = std::max(std::max(a,b),c);
}
template <int dim,int spacedim, int chartdim>
Point<chartdim> Ellipsoid<dim,spacedim,chartdim>::pull_back(const Point<spacedim> &space_point) const
{
double x,y,z, u,v,w;
// get point on ellipsoid
x = space_point[0];
y = space_point[1];
z = space_point[2];
std::cout << "using a,b,c: " << std::endl;
std::cout << a << " " << b << " " << c << std::endl;
std::cout << "from pull_back: " << std::endl;
std::cout << "space_point: " << std::endl;
std::cout << x << " " << y << " " << z << std::endl;
// map ellipsoid point onto sphere
u = x/a;
v = y/b;
w = z/c;
std::cout << "pulls back to : " << std::endl;
std::cout << u << " " << v << " " << w << std::endl;
std::cout << "on sphere." << std::endl;
Point<spacedim> p(u,v,w);
// use reference_sphere's pull_back function
Point<spacedim> q = reference_sphere.pull_back(p);
Point<chartdim> chart_point;
std::cout << "sphere pull_back: " << std::endl;
std::cout << q[0] << " " << q[1] << " " << q[2] << std::endl;
std::cout << "r theta phi" << std::endl;
std::cout << "..........." << std::endl;
chart_point[0] = q[1];
chart_point[1] = q[2];
// return (theta,phi) in the chart domain
return chart_point;
}
template <int dim,int spacedim, int chartdim>
Point<spacedim> Ellipsoid<dim,spacedim,chartdim>::push_forward(const Point<chartdim> &chart_point) const
{
double theta,phi, x,y,z;
phi = chart_point[0];
theta = chart_point[1];
Point<spacedim> p(max_axis,theta,phi);
// map theta,phi in chart domain onto reference_sphere with radius max_axis
Point<spacedim> X = reference_sphere.push_forward(p);
// map point on sphere onto ellipsoid
x = a*X[0];
y = b*X[1];
z = c*X[2];
Point<spacedim> space_point(x,y,z);
// return point on ellipsoid
return space_point;
}
Point<3> grid_transform (const Point<3> &X)
{
// transform points on sphere onto ellipsoid
double a,b,c;
a = 1.0; b = 1.0; c = 1.0;
double x,y,z;
x = a*X(0);
y = b*X(1);
z = c*X(2);
return Point<3>(x,y,z);
}
void assemble_mesh_and_manifold()
{
const int dim = 2;
const int spacedim = 3;
const int chartdim = 2;
double a,b,c;
a = 1; b=1; c=1;
Ellipsoid<dim,spacedim,chartdim> ellipsoid(a,b,c);
Triangulation<dim,spacedim> tria;
// generate coarse spherical mesh
GridGenerator::hyper_sphere (tria, Point<spacedim>(0.0,0.0,0.0), 1.0);
for (Triangulation<dim,spacedim>::active_cell_iterator cell=tria.begin_active(); cell!=tria.end(); ++cell)
cell->set_all_manifold_ids(0);
print_mesh_info(tria, "spherical_mesh.vtk");
GridTools::transform(&grid_transform, tria);
tria.set_manifold(0,ellipsoid);
tria.refine_global(1);
print_mesh_info(tria, "ellipsoidal_mesh.vtk");
// transform the mesh to be ellipsoidal
//GridTools::transform(std_cxx11::bind(&Ellipsoid<dim,spacedim>::push_forward,std_cxx11::cref(ellipsoid),std_cxx11::_1), tria);
//print_mesh_info(tria, "ellipsoidal_mesh.vtk");
//GridOut grid_out;
//std::ofstream out ("ellipsoidal_mesh.vtk");
//grid_out.write_vtk (tria, out);
////std::cout << " written to " << filename
// << std::endl
}
// @sect3{The main function}
// Finally, the main function. There isn't much to do here, only to call the
// subfunctions.
int main ()
{
assemble_mesh_and_manifold();
}
ellipsoidal_mesh.vtk
Description: Binary data
