Dear all,
I'm am trying to reproduce with my implementation, the results in the
Photonic Crystal computations performed in [1]. Here, the author uses a
grid with an inner disk with radius R=0.475, and for FEM it is used the
software Concepts [2] that implements curvilinear elements denoted Blending
technique, or transfinite interpolation in quadrilaterals [3], [4,p. 144].
Exponential convergence was obtained for the p-FEM version.
By using the deal.II SphericalManifold, I can only get with the attached
code up to R<0.46, before I get the exception:
"The image of the mapping applied to cell with center [0.446317 -0.206317]
is distorted. The cell geometry or the mapping are invalid, giving a
non-positive volume fraction of -0.000129188 in quadrature point 54."
I am attaching a minimal test case (modified from step-10) that reproduces
the errors:
mesh_test.cc
and the code for the same mesh used in [1], adapted to deal.II standards in:
unit_cell.h
The question is then if it is possible to get this example running for
R=0.475, in order to reproduce the results in [1]. Is there something I am
doing wrong that can be improved? Can this behaviour be explained by round
off's?
Furthermore, is anyone implementing transfinite interpolation [3,4] in
deal.II?
Thanks in advance, any help is greatly appreciated!
References:
[1] Engström, C.,Spectral approximation of quadratic operator polynomials
arising in photonic band structure calculations, Numer. Math. (2014) 126:
413.
[2] Concepts web page:
http://wiki.math.ethz.ch/Concepts/WebRoot?redirectedfrom=Concepts.WebHome
[3] W. J. Gordon, C. A. Hall, Construction of curvilinear coordinate
systems and application to mesh generation Int. J. Num. Meth. Eng., Vol. 7
(1973), pp. 461-477
[4] Pavel Solin, Karel Segeth, Ivo Dolezel, Higher-Order Finite Element
Methods, Studies in Advanced Mathematics, CRC Press, 2003.
--
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 [email protected].
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
*/
#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>
// added lines
#include "unit_cell.h"
#include <deal.II/numerics/data_out.h>
namespace Step10
{
using namespace dealii;
const long double pi = 3.141592653589793238462643;
template <int dim>
void gnuplot_output()
{
std::cout << "Output of grids into gnuplot files:" << std::endl
<< "===================================" << std::endl;
Triangulation<dim> triangulation;
GridGenerator::hyper_ball (triangulation);
static const SphericalManifold<dim> boundary;
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold (0, boundary);
for (unsigned int refinement=0; refinement<2;
++refinement, triangulation.refine_global(1))
{
std::cout << "Refinement level: " << refinement << std::endl;
std::string filename_base = "ball";
filename_base += '0'+refinement;
for (unsigned int degree=1; degree<4; ++degree)
{
std::cout << "Degree = " << degree << std::endl;
const MappingQ<dim> mapping (degree);
GridOut grid_out;
GridOutFlags::Gnuplot gnuplot_flags(false, 30);
grid_out.set_flags(gnuplot_flags);
std::string filename = filename_base+"_mapping_q";
filename += ('0'+degree);
filename += ".dat";
std::ofstream gnuplot_file (filename.c_str());
grid_out.write_gnuplot (triangulation, gnuplot_file, &mapping);
}
std::cout << std::endl;
}
}
template <int dim>
void compute_pi_by_area ()
{
std::cout << "Computation of Pi by the area:" << std::endl
<< "==============================" << std::endl;
// const QGauss<dim> quadrature(4);
for (double R=0.45; R < 0.5; R += 0.005)
for (unsigned int degree=1; degree<11; ++degree)
{
std::cout << "Radius = " << R << std::endl;
std::cout << "Degree = " << degree << std::endl;
Triangulation<dim> triangulation;
// GridGenerator::hyper_ball (triangulation);
// new lines
Geom_parameters gp;
unsigned int curved_faces_label=10;
cell_rod ( triangulation, 1.0, R, gp, false );
static const SphericalManifold<dim> boundary;
triangulation.set_manifold ( curved_faces_label, boundary ); // labels defined with the grid
// triangulation.set_all_manifold_ids_on_boundary (0);
const MappingQ<dim> mapping (degree, true);
const QGauss<dim> quadrature(degree);
// const FE_Q<dim> dummy_fe (1);
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;
// No refinements
for (unsigned int refinement=0; refinement<1;
++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) {
if (r < R)
area += fe_values.JxW (i);
}
};
double new_area = area/(R*R);
table.add_value("eval.pi", static_cast<double> (new_area));
table.add_value("error", static_cast<double> (std::fabs(new_area-pi)));
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
Vector<double> dummy ( dof_handler.n_dofs() );
string name = "grid";
data_out.add_data_vector (dummy, name );
std::ostringstream ss;
ss <<"grid_R" << R << "_p" << degree << "_ref" << refinement << ".vtk";
std::ofstream output (ss.str().c_str());
data_out.build_patches (mapping, 8, 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;
}
#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/tria_boundary_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_reordering.h>
// #include <deal.II/grid/manifold_lib.h>
#include <fstream>
#include <assert.h>
#include <map>
#define PI 3.14159265358979323846
using namespace dealii;
using namespace std;
struct Geom_parameters {
/*
* model: geometric configuration, parameters must be dependent on the model problem
*/
unsigned int n_scatterers, scatterers_coarse_cells, n_balls, coated_coarse_cells, defective_index;
double coat_radius;
bool coated=false;
bool defect=false;
bool defect_n2=false;
vector< Point<2> > ball_centers;
vector<double> radius;
vector<double> radius_PML;
vector<unsigned int> ball_labels;
// double a;
string name;
};
// In the article "Spectral approximation of quadratic operator polynomials arising in photonic band structure calculations", Christian Engström 2012,
//
// This mesh is used with R=0.475 and exponential convergence is obtained.
// Check figure 1, and figure 3 on the mentioned paper.
//
// http://link.springer.com/article/10.1007/s00211-013-0568-y
//
// Software Concepts:
// http://wiki.math.ethz.ch/Concepts/WebRoot?redirectedfrom=Concepts.WebHome
void cell_rod ( Triangulation<2> &tria, double L, double r, // 2*r/L < 1.0
Geom_parameters &gp, bool mesh_out)
{
double l = 0.5*L, d = 0.5*r, q=1.0/sqrt(2.0); // q: corner points factor
const unsigned int n = 2, n_vert = 24+1, n_cell = 19+1;
bool info = false, plot_eps = true; // for printing cells, faces and vertex values
int fv = 4, fc = 1; // f: fixed, displacement of index on the vertexes and cells ... inner objects
// geometric information-----------------------
gp.n_scatterers = 1;
gp.scatterers_coarse_cells = 12;
gp.n_balls = 1; // just one Manifold ... everything centered at origin
gp.ball_centers.resize( gp.n_balls );
gp.radius.resize(gp.n_balls);
gp.radius[0]=r;
gp.ball_centers[0] = Point<2>( 0.0, 0.0 );
gp.name = "cell_rod";
if(mesh_out) printf("%% defined by using %d balls\n", gp.n_balls);
std::vector<Point<2> > vertices(n_vert);
vertices[0] = Point<2>( 0.0, 0.0 );
// inner 4 squares
vertices[1] = Point<2>( d, 0.0 );
vertices[2] = Point<2>( .5*d, .5*d );
vertices[3] = Point<2>( 0.0, d );
vertices[4] = Point<2>( -.5*d, .5*d );
vertices[5] = Point<2>( -d, 0.0 );
vertices[6] = Point<2>( -.5*d, -.5*d );
vertices[7] = Point<2>( 0.0, -d );
vertices[8] = Point<2>( .5*d, -.5*d );
// touching circle
vertices[9] = Point<2>( r, 0.0 );
vertices[10] = Point<2>( r*q, r*q );
vertices[11] = Point<2>( 0.0, r );
vertices[12] = Point<2>( -r*q, r*q );
vertices[13] = Point<2>( -r, 0.0 );
vertices[14] = Point<2>( -r*q, -r*q );
vertices[15] = Point<2>( 0.0, -r );
vertices[16] = Point<2>( r*q, -r*q );
vertices[17] = Point<2>( l, 0.0 );
vertices[18] = Point<2>( l, l );
vertices[19] = Point<2>( 0.0, l );
vertices[20] = Point<2>( -l, l );
vertices[21] = Point<2>( -l, 0.0 );
vertices[22] = Point<2>( -l, -l );
vertices[23] = Point<2>( 0.0, -l );
vertices[24] = Point<2>( l, -l );
//vertices[] = Point<2>( , );
std::vector< std::vector<int> > cell_v (n_cell, std::vector<int>(4));
cell_v[0][0] = 0; // cell, i = 0
cell_v[0][1] = 8;
cell_v[0][2] = 2;
cell_v[0][3] = 1;
cell_v[1][0] = 0; // cell, i = 1
cell_v[1][1] = 2;
cell_v[1][2] = 4;
cell_v[1][3] = 3;
cell_v[2][0] = 0; // cell, i = 2
cell_v[2][1] = 4;
cell_v[2][2] = 6;
cell_v[2][3] = 5;
cell_v[3][0] = 0; // cell, i = 3
cell_v[3][1] = 6;
cell_v[3][2] = 8;
cell_v[3][3] = 7;
cell_v[4][0] = 8; // cell, i = 4
cell_v[4][1] = 16;
cell_v[4][2] = 1;
cell_v[4][3] = 9;
cell_v[5][0] = 2; // cell, i = 5
cell_v[5][1] = 1;
cell_v[5][2] = 10;
cell_v[5][3] = 9;
cell_v[6][0] = 2; // cell, i = 6
cell_v[6][1] = 10;
cell_v[6][2] = 3;
cell_v[6][3] = 11;
cell_v[7][0] = 4; // cell, i = 7
cell_v[7][1] = 3;
cell_v[7][2] = 12;
cell_v[7][3] = 11;
cell_v[8][0] = 4; // cell, i = 8
cell_v[8][1] = 12;
cell_v[8][2] = 5;
cell_v[8][3] = 13;
cell_v[9][0] = 6; // cell, i = 9
cell_v[9][1] = 5;
cell_v[9][2] = 14;
cell_v[9][3] = 13;
cell_v[10][0] = 6; // cell, i = 10
cell_v[10][1] = 14;
cell_v[10][2] = 7;
cell_v[10][3] = 15;
cell_v[11][0] = 8; // cell, i = 11
cell_v[11][1] = 7;
cell_v[11][2] = 16;
cell_v[11][3] = 15;
// cell 12
cell_v[12][0] = 16; // cell, i = 12
cell_v[12][1] = 24;
cell_v[12][2] = 9;
cell_v[12][3] = 17;
// cell 13
cell_v[13][0] = 10; // cell, i = 13
cell_v[13][1] = 9;
cell_v[13][2] = 18;
cell_v[13][3] = 17;
// cell 14
cell_v[14][0] = 10; // cell, i = 14
cell_v[14][1] = 18;
cell_v[14][2] = 11;
cell_v[14][3] = 19;
// cell 15
cell_v[15][0] = 12; // cell, i = 15
cell_v[15][1] = 11;
cell_v[15][2] = 20;
cell_v[15][3] = 19;
// cell 16
cell_v[16][0] = 12; // cell, i = 16
cell_v[16][1] = 20;
cell_v[16][2] = 13;
cell_v[16][3] = 21;
// cell 17
cell_v[17][0] = 14; // cell, i = 17
cell_v[17][1] = 13;
cell_v[17][2] = 22;
cell_v[17][3] = 21;
// cell 18
cell_v[18][0] = 14; // cell, i = 18
cell_v[18][1] = 22;
cell_v[18][2] = 15;
cell_v[18][3] = 23;
// cell 19
cell_v[19][0] = 16; // cell, i =
cell_v[19][1] = 15;
cell_v[19][2] = 24;
cell_v[19][3] = 23;
if (info) {
for (unsigned int i = 0; i < n_cell; ++i)
if(info) printf("cell[%d] = (%d,%d,%d,%d)\n", i, cell_v[i][0], cell_v[i][1],cell_v[i][2],cell_v[i][3]);
}
std::vector<CellData<2> > cells (n_cell, CellData<2>());
unsigned int cell_idx = 0;
for (unsigned int i=0; i < n_cell; ++i) {
// printf("%d: ", i);
for (unsigned int j=0; j < 4; ++j) {
cells[i].vertices[j] = cell_v[i][j];
if(info) printf("%d ",cell_v[i][j]);
}
if(info) printf(";\n");
cell_idx++;
}
for (unsigned int i=0; i < n_vert; ++i) {
if(info) printf(" %f, %f;\n", vertices[i][0],vertices[i][1] );
}
if(info) printf("\n\n");
//printf("after grid ordering\n");
tria.create_triangulation (
std::vector<Point<2> >(&vertices[0], &vertices[n_vert]),
cells,
SubCellData()); // no boundary information
if(info) printf("after creating triangulation\n");
// colorizing ------------------------------------------------------------------------------------
double eps = 1e-5, mid;
unsigned int label=10;
// cell way
for (Triangulation<2>::active_cell_iterator
cell = tria.begin_active(); cell!=tria.end(); ++cell)
{
// cell->set_all_manifold_ids (label);
for (unsigned int f=0;
f < GeometryInfo<2>::faces_per_cell;++f) {
const Point<2> face_center = cell->face(f)->center(),
p0 = cell->face(f)->vertex(0) , p1 = cell->face(f)->vertex(1);
double d0 = p0.distance( Point<2>( 0.0,0.0 ) ), d1 = p1.distance( Point<2>( 0.0,0.0 ) );
if(info) printf(" f=%d, dist=%f\n", f, cell->face(f)->center().distance( Point<2>( 0.0,0.0 ) ) );
if( (cell->face(f)->center().distance( Point<2>( 0.0,0.0 )) - mid) < eps ) {
if(info) printf(" ... set one! ... dist=%f\n",
cell->face(f)->center().distance( Point<2>( 0.0,0.0 ) ) );
cell->face(f)->set_all_manifold_ids (label+1);
} // if face of circle
if ( abs((d0+d1)-2*r) < eps ) { // disk faces
cell->face(f)->set_all_manifold_ids (label);
} // outer faces
} // face
if ( cell->center().distance( Point<2>( 0.0,0.0 ) ) < eps ) { // eliminating origin
cell->set_all_manifold_ids (label+1);
}
} // cell
// end-colorizing --------------------------------------------------------------------------------
if(mesh_out) { // printing the resulting geometry
//printf("printing original mesh to eps ... \n");
std::ofstream gridfile ("coarse.eps");
GridOut grid_out;
grid_out.write_eps (tria, gridfile);
// printing vertices
printf("vert=[\n");
for (unsigned int j = 0; j < n_vert; j++) {
printf("\t%.3f, %.3f;\n", vertices[j](0), vertices[j](1));
}
printf("];\n\n");
// printing cell indexes
printf("cells=[\n");
for (unsigned int j = 0; j < n_cell; j++) {
printf("\t%d, %d,%d, %d;\n", cell_v[j][0], cell_v[j][1], cell_v[j][2], cell_v[j][3] );
}
printf("];\n");
}
}
void cell_rod ( Triangulation<2> &tria, double L, double diameter_ratio) {
Geom_parameters gp;
cell_rod ( tria, L, diameter_ratio, gp, false);
}
void cell_rod ( Triangulation<2> &tria, double L, double diameter_ratio,
Geom_parameters &gp) {
cell_rod ( tria, L, diameter_ratio, gp, false);
}
