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, 5/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;
// Description of the upper face of the square
class UpperManifold: public ChartManifold<2,2,1> {
public:
virtual Point<1>
pull_back(const Point<2> &space_point) const;
virtual Point<2>
push_forward(const Point<1> &chart_point) const;
};
Point<1> UpperManifold::pull_back(const Point<2> &space_point) const {
Point<1> p(space_point[0]);
return p;
}
Point<2> UpperManifold::push_forward(const Point<1> &chart_point) const {
double x = chart_point[0];
// return Point<2> (x, -2*x*x + 2*x + 1); // Parabole
return Point<2> ( x, 0.25*sin ( PI*x) + 1); // Trigonometric
}
template <int dim>
void compute_pi_by_area ()
{
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;
Triangulation<dim> triangulation;
// new lines
unsigned int curved_faces_label=3; // see GridGenerator::hyper_rectangle, colorize=true
GridGenerator::hyper_cube (triangulation, 0.0, 1.0, true);
static const UpperManifold boundary;
triangulation.set_manifold ( curved_faces_label, boundary );
const MappingQ<dim> mapping (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<5;
++refinement, triangulation.refine_global (1))
{
table.add_value("cells", triangulation.n_active_cells());
dof_handler.distribute_dofs (dummy_fe);
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 new_area = 0.5/(area-1.0); // not working
table.add_value("eval.pi", static_cast<double> (new_area));
table.add_value("error", static_cast<double> (std::fabs(new_area-PI)));
if (refinement == 0) {
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);
}
};
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;
};
}
}
int main ()
{
try
{
std::cout.precision (16);
Step10::compute_pi_by_area<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;
}
