In the same spirit I would like to share a code that uses ChartManifold,
for generating a mapping of an ellipse applied to HyperBall.
El miércoles, 5 de abril de 2017, 14:47:17 (UTC+2), Juan Carlos Araujo
Cabarcas escribió:
>
> Dear all,
> I recently found the thread: Something wrong with ChartManifold when
> chartdim=1:
> https://groups.google.com/forum/#!topic/dealii/Sfo9xKoeRpw
> a more straight answer to this thread.
>
> I took the liberty to copy David's code and modify it, so that the upper
> face of a square follows: y(x)=0.25*sin(PI*x) + 1.
> The function is not the same as I first proposed, but from this example it
> is pretty clear how to proceed.
> The result is plotted in the attached .png, along with a variation from
> step-10 in the deal.ii tutorial.
>
> I hope this may be of any use to somebody!
>
> Juan Carlos Araújo,
> Umeå Universitet
> El martes, 9 de febrero de 2016, 17:10:51 (UTC+1), Juan Carlos Araujo
> Cabarcas escribió:
>>
>> Dear all,
>>
>> I am working with the wave equation with variable refractive indexes
>> within the domain.
>> I receive meshes that come from a shape optimization routine, and my task
>> is to
>> run on arbitrarily curved elements resulting from the optimization.
>> Step-10 and others show how to use mappings and Manifolds to curve
>> elements following basic shapes like circles, cylinders etc.
>>
>> My guess is that it should be possible in deal.II to define my own
>> Manifolds based on how I want to bend my edges (optimization routine).
>>
>> I have prepared a very basic example on the matlab code attached that has
>> the info to plot the attached figure. There, the edges follow the equation
>> |x|^p + |y|^p=1, with p=6.
>>
>> It would be very illustrative to apply a mapping like what is done in
>> step-10 on a circle, but for the presented case. I would appreciate any
>> advise on how to achieve this, hopefully wth this example =)
>>
>> Thanks in advance!
>>
>> Juan Carlos Araújo-Cabarcas
>> Umeå Universitet.
>>
>>
>>
--
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.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2001 - 2014 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.
*
* ---------------------------------------------------------------------
*
* Authors: Wolfgang Bangerth, Ralf Hartmann, University of Heidelberg, 2001
* Modified from step-10 by Juan Carlos Araújo, 7/apr/2017
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>
#include <iostream>
#include <fstream>
#include <cmath>
#include <deal.II/numerics/data_out.h>
#define PI 3.14159265358979323846
namespace Step10
{
using namespace dealii;
/* ************** READ ME ****************************************************
Illustration of the use of ChartManifold, for generating a mapping
of an ellipse starting from HyperBall.
Juan Carlos Araújo, 7/apr/2017
*/
Point<2> to_manifold (Point<2> space_point, double a, double b) {
Point<2> R = Point<2>(space_point[0],space_point[1]), p;
double rho, theta;
theta = atan2(R[1],R[0]); // theta
if (theta < 0)
theta += 2*PI;
p = Point<2> ( a*cos(theta), b*sin(theta) );
return p;
}
// Rewritten from SphericalManifold 2D
class EllipseManifold : public ChartManifold<2, 2, 2> {
public:
EllipseManifold(double a, double b);
virtual Point<2>
pull_back(const Point<2> &space_point) const;
virtual Point<2>
push_forward(const Point<2> &chart_point) const;
virtual Point<2>
get_new_point(const Quadrature<2> &quad) const;
double a, b;
private:
static Point<2> get_periodicity();
};
EllipseManifold::EllipseManifold(double a, double b):
ChartManifold<2,2,2>(EllipseManifold::get_periodicity()), a(a),b(b)
{}
Point<2>
EllipseManifold::get_periodicity()
{
Point<2> periodicity;
periodicity[2-1] = 2*PI;
return periodicity;
}
Point<2>
EllipseManifold::get_new_point(const Quadrature<2> &quad) const {
return ChartManifold<2,2,2>::get_new_point(quad);
}
Point<2>
EllipseManifold::push_forward(const Point<2> &spherical_point) const {
Assert(spherical_point[0] >=0.0,
ExcMessage("Negative radius for given point."));
const double rho = spherical_point[0];
const double theta = spherical_point[1];
Point<2> p;
if (rho > 1e-10)
p[0] = a*cos(theta);
p[1] = b*sin(theta);
return p;
}
Point<2>
EllipseManifold::pull_back(const Point<2> &space_point) const {
Point<2> R = Point<2>(space_point[0],space_point[1]), p;
double rho;
p[1] = atan2(R[1]/b,R[0]/a); // theta in the frame of the Ellipse
if (p[1] < 0)
p[1] += 2*PI;
R = Point<2> ( a*cos(p[1]), b*sin(p[1]) );
rho = R.norm();
p[0] = rho;
return p;
}
template <int dim>
void create_grid () {
std::cout << "Computation of Pi:" << std::endl
<< "==============================" << std::endl;
// The case p=1, is trivial: no bending!
for (unsigned int degree=2; degree<11; ++degree) {
std::cout << "Degree = " << degree << std::endl;
double axisX=1.5, axisY=1.0;
Triangulation<dim> triangulation;
static const EllipseManifold boundary (axisX, axisY);
GridGenerator::hyper_ball (triangulation);
triangulation.set_manifold ( 0, boundary );
// Here, we move the outer nodes of HyperBall to coincide with the Ellipse
Triangulation<2>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell) {
for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i) {
Point<2> &v = cell->vertex(i), tp;
if ( abs (v.square () - 1.0) < 1e-10 ) {
Point<2> t = to_manifold (v,axisX, axisY);
v(0) = t(0); v(1) = t(1);
}
}
}
const MappingQ<dim> mapping ( degree, true ); // degree, true
const QGauss<dim> quadrature(degree);
const FE_Q<dim> dummy_fe (QGaussLobatto<1>(degree + 1));
DoFHandler<dim> dof_handler (triangulation);
FEValues<dim> fe_values (mapping, dummy_fe, quadrature,
update_JxW_values);
ConvergenceTable table;
for (unsigned int refinement=0; refinement<4;
++refinement, triangulation.refine_global (1))
{
table.add_value("cells", triangulation.n_active_cells());
dof_handler.distribute_dofs (dummy_fe);
// outputting grid before using FE ...
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
Vector<double> dummy ( dof_handler.n_dofs() );
std::string name = "grid";
data_out.add_data_vector (dummy, name );
std::ostringstream ss;
ss <<"grid_p" << degree << "_ref" << refinement << ".vtk";
std::ofstream output (ss.str().c_str());
data_out.build_patches (mapping, 20, DataOut<dim>::curved_inner_cells);
data_out.write_vtk (output);
long double area = 0;
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell) {
double r = cell->center ().distance(Point<dim>(0,0));
fe_values.reinit (cell);
for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i) {
area += fe_values.JxW (i);
}
};
double pi_approx = area/(axisX*axisY);
table.add_value("eval.pi", static_cast<double> (pi_approx));
table.add_value("error", static_cast<double> (std::fabs(pi_approx-PI)));
table.omit_column_from_convergence_rate_evaluation("cells");
table.omit_column_from_convergence_rate_evaluation("eval.pi");
table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
table.set_precision("eval.pi", 16);
table.set_scientific("error", true);
table.write_text(std::cout);
std::cout << std::endl;
} // ref
} // degree
} // create_grid
}
int main ()
{
try
{
std::cout.precision (16);
Step10::create_grid <2> ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
