Dear all,
Since the inclusion of Transfinite interpolation, I have been successful on
working with this powerful technique in my research. I had coded a mesh
implementing concentric circles, where the inner most is shifted a small
distance s. All concentric circles are labeled 100+i, where i is in a loop
on all circles.
The mesh was working after 4/Apr/17 when I added the following entry in the
forum:
https://groups.google.com/forum/#!topic/dealii/hCZqv9g6mKk
I installed the development version of dealii on the 17/nov/17. The details
of the installation are also in the forum:
https://groups.google.com/forum/#!topic/dealii/ee2w2X987ig
and recently I went back and tried to run thhis code, obtaining the error
at the end of the email.
I prepared a minimal example that includes the definition of my grid. It
includes how I colorize each face as explained before, and the rest are
colorized "1" as the documentation suggests.
The symptoms are the following:
- If I use refinement=0, everything works and the mesh looks quite nice in
zeros.vtk and surface with lines mode in Paraview.
- For any higher refinement, it gives the error below.
Maybe there has been a major change in TransfiniteInterpolation since
Apr/2017, or maybe the way I colorize is not good anymore for some weird
reason. Any ideas or comments would be greatly appreciated.
Juan Carlos Araújo,
Umeå Universitet.
terminate called after throwing an instance of 'dealii::Mapping<2,
2>::ExcTransformationFailed'
what():
--------------------------------------------------------
An error occurred in line <1554> of file
</path/dealii/source/grid/manifold_lib.cc> in function
typename dealii::Triangulation<dim, spacedim>::cell_iterator
dealii::TransfiniteInterpolationManifold<dim,
spacedim>::compute_chart_points(const dealii::ArrayView<const
dealii::Point<spacedim> >&, dealii::ArrayView<dealii::Point<dim> >) const
[with int dim = 2; int spacedim = 2; typename dealii::Triangulation<dim,
spacedim>::cell_iterator = dealii::TriaIterator<dealii::CellAccessor<2, 2>
>]
The violated condition was:
false
Additional information:
--------------------------------------------------------
Aborted (core dumped)
--
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) 1999 - 2013 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.
*
* ---------------------------------------------------------------------
*
*/
#include <deal.II/grid/tria.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/tria_boundary_lib.h>
#include <deal.II/grid/manifold_lib.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/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/lac/constraint_matrix.h>
// #include "../../grid/arbitrary_resonators.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <typeinfo>
#include <limits>
#include <iomanip>
using namespace dealii;
using namespace std;
//---------------------------- GRID DEFINITION ---------------------------------
struct Geom_parameters {
unsigned int n_scatterers, scatterers_coarse_cells,
n_balls, coated_coarse_cells, defective_index;
unsigned int manifold_label,layers_label,bc_label;
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;
string name;
};
void concentric_disks ( Triangulation<2> &tria, double s, vector<double> x,
Geom_parameters &gp, bool mesh_out)
{
double r = x[0], d = 0.5*x[0], q=1.0/sqrt(2.0); //l = x[1], q: corner points factor
const unsigned int nlay = x.size()-1, vlayer=8, lcells=8,
n_vert = (2+nlay)*vlayer+1, n_cell = 3+(1+nlay)*lcells+1; // 19
bool info = false; //, plot_eps = true, for printing cells, faces and vertex values
int fc = 4; //fv = 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 = x.size(); // the resonator and all other centered at the origin
gp.ball_centers.resize( gp.n_balls );
gp.layers_label=1;
gp.manifold_label=10;
gp.bc_label=0;
gp.radius.resize(gp.n_balls);
gp.radius[0]=x[0];
gp.ball_centers[0] = Point<2>( s, 0.0 );
for (unsigned int k = 1; k < x.size(); k++) {
gp.radius[k]=x[k];
gp.ball_centers[k] = Point<2>( 0.0, 0.0 );
}
gp.name = "concentric_disks";
if(mesh_out) printf("%% defined by using %d balls\n", gp.n_balls);
if (info) printf (" Entering %s ... \n", gp.name.c_str () );
std::vector<Point<2> > vertices(n_vert);
vertices[0] = Point<2>( 0.0+s, 0.0 );
// inner 4 squares
vertices[1] = Point<2>( d+s, 0.0 );
vertices[2] = Point<2>( .5*d+s, .5*d );
vertices[3] = Point<2>( 0.0+s, d );
vertices[4] = Point<2>( -.5*d+s, .5*d );
vertices[5] = Point<2>( -d+s, 0.0 );
vertices[6] = Point<2>( -.5*d+s, -.5*d );
vertices[7] = Point<2>( 0.0+s, -d );
vertices[8] = Point<2>( .5*d+s, -.5*d );
// touching circle
for (unsigned int k=0; k < x.size();k++) {
double z = (k==0)?s:0.0;
r = x[k];
vertices[8+k*vlayer+1] = Point<2>( r+z, 0.0 );
vertices[8+k*vlayer+2] = Point<2>( r*q+z, r*q );
vertices[8+k*vlayer+3] = Point<2>( 0.0+z, r );
vertices[8+k*vlayer+4] = Point<2>( -r*q+z, r*q );
vertices[8+k*vlayer+5] = Point<2>( -r+z, 0.0 );
vertices[8+k*vlayer+6] = Point<2>( -r*q+z, -r*q );
vertices[8+k*vlayer+7] = Point<2>( 0.0+z, -r );
vertices[8+k*vlayer+8] = Point<2>( r*q+z, -r*q );
}
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;
unsigned int m;
// layer cells
for (unsigned int k = 1; k < x.size(); k++) {
m=k+1;
// cell 12
cell_v[fc+k*lcells+0][0] = 7+8*k+1; //16; // cell, i = 12
cell_v[fc+k*lcells+0][1] = 7+8*m+1; //24;
cell_v[fc+k*lcells+0][2] = 0+8*k+1; // 9;
cell_v[fc+k*lcells+0][3] = 0+8*m+1; //17;
// cell 13
cell_v[fc+k*lcells+1][0] = 1+8*k+1; //10; // cell, i = 13
cell_v[fc+k*lcells+1][1] = 0+8*k+1; // 9;
cell_v[fc+k*lcells+1][2] = 1+8*m+1; //18;
cell_v[fc+k*lcells+1][3] = 0+8*m+1; //17;
// cell 14
cell_v[fc+k*lcells+2][0] = 1+8*k+1; //10; // cell, i = 14
cell_v[fc+k*lcells+2][1] = 1+8*m+1; //18;
cell_v[fc+k*lcells+2][2] = 2+8*k+1; //11;
cell_v[fc+k*lcells+2][3] = 2+8*m+1; //19;
// cell 15
cell_v[fc+k*lcells+3][0] = 3+8*k+1; //12; // cell, i = 15
cell_v[fc+k*lcells+3][1] = 2+8*k+1; //11;
cell_v[fc+k*lcells+3][2] = 3+8*m+1; //20;
cell_v[fc+k*lcells+3][3] = 2+8*m+1; //19;
// cell 16
cell_v[fc+k*lcells+4][0] = 3+8*k+1; //12; // cell, i = 16
cell_v[fc+k*lcells+4][1] = 3+8*m+1; //20;
cell_v[fc+k*lcells+4][2] = 4+8*k+1; //13;
cell_v[fc+k*lcells+4][3] = 4+8*m+1; //21;
// cell 17
cell_v[fc+k*lcells+5][0] = 5+8*k+1; //14; // cell, i = 17
cell_v[fc+k*lcells+5][1] = 4+8*k+1; //13;
cell_v[fc+k*lcells+5][2] = 5+8*m+1; //22;
cell_v[fc+k*lcells+5][3] = 4+8*m+1; //21;
// cell 18
cell_v[fc+k*lcells+6][0] = 5+8*k+1; //14; // cell, i = 18
cell_v[fc+k*lcells+6][1] = 5+8*m+1; //22;
cell_v[fc+k*lcells+6][2] = 6+8*k+1; //15;
cell_v[fc+k*lcells+6][3] = 6+8*m+1; //23;
// cell 19
cell_v[fc+k*lcells+7][0] = 7+8*k+1; //16; // cell, i =
cell_v[fc+k*lcells+7][1] = 6+8*k+1; //15;
cell_v[fc+k*lcells+7][2] = 7+8*m+1; //24;
cell_v[fc+k*lcells+7][3] = 6+8*m+1; //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");
double eps = 1e-5*x[0];
unsigned int label=100; // layers=1, origin_cells=100,
for (Triangulation<2>::active_cell_iterator
cell = tria.begin_active(); cell!=tria.end(); ++cell)
{
cell->set_all_manifold_ids (1); // not faces ... for Transfinite interpolation
for (unsigned int k=0; k < gp.n_balls; ++k) {
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( gp.ball_centers[k] ), d1 = p1.distance( gp.ball_centers[k] );
if(info) printf(" f=%d, dist=%f, \t r=%f\n", f,
cell->face(f)->center().distance( gp.ball_centers[k] ),gp.radius[k] );
if( (abs( d0 - gp.radius[k]) < eps ) &&
(abs( d1 - gp.radius[k]) < eps )) {
cell->face(f)->set_all_manifold_ids (label+k);
if(info) printf("Manifold of ball %d! dist=%f, \t r=%f\n", k,
cell->face(f)->center().distance( gp.ball_centers[k] ), gp.radius[k] );
}
}
} // face
} // cell
// end-colorizing --------------------------------------------------------------------------------
// print grid
if(mesh_out) {
std::ofstream gridfile ("grid_cc.eps");
GridOut grid_out;
grid_out.write_eps (tria, gridfile);
std::ostringstream ss; // auxiliar string used for naming files
ss << "coarse_grid_" << "cc" << ".m";
std::ofstream output (ss.str().c_str());
output << "vert=[" << endl;
for (unsigned int j = 0; j < n_vert; j++) {
//printf("\t%.3f, %.3f;\n", vertices[j](0), vertices[j](1));
output << "\t" << vertices[j](0) << ",\t" << vertices[j](1) << ";" << endl;
output.flush();
}
output << "\n];" << endl;
output << "cells=[" << endl;
for (unsigned int j = 0; j < n_cell; j++) {
output << "\t" << cell_v[j][0] << ",\t" << cell_v[j][1]
<< ",\t" << cell_v[j][2] << ",\t" << cell_v[j][3] <<";" << endl;
output.flush();
}
output << "\n];" << endl;
output.close();
}
if (info) printf (" Exit %s ... \n", gp.name.c_str () );
}
void concentric_disks ( Triangulation<2> &tria, vector<double> x) {
Geom_parameters gp;
concentric_disks ( tria, 0.0, x, gp, false);
}
void concentric_disks ( Triangulation<2> &tria, vector<double> x,
Geom_parameters &gp) {
concentric_disks ( tria, 0.0, x, gp, false);
}
//--------------------- MAIN CLASS ---------------------------------------------
template <int dim>
class Mygrid
{
public:
Mygrid (unsigned int p, unsigned int r);
void run ();
private:
void make_grid ();
void setup_system();
Geom_parameters gp;
vector< PolarManifold<dim> > balls;
TransfiniteInterpolationManifold<dim> inner_manifold;
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
MappingQGeneric<dim> mapping;
ConstraintMatrix constraints;
unsigned int degree, refinement;
};
template <int dim>
Mygrid<dim>::Mygrid (unsigned int p, unsigned int r)
:
fe (QGaussLobatto<1>(p + 1)),
dof_handler (triangulation),
mapping (p),
degree (p), refinement (r)
{
printf ("Degree=%d, refinement=%d\n", degree, refinement);
}
template <int dim>
void Mygrid<dim>::make_grid () {
double s=0.1;
std::vector<double> x;
x.push_back(1.0);
x.push_back(1.5);
x.push_back(2.0);
x.push_back(2.5);
x.push_back(3.0);
concentric_disks ( triangulation, s, x, gp, true);
balls.resize(gp.n_balls);
for (unsigned int i = 0; i < gp.n_balls; i++) {
new (&balls[i]) PolarManifold<2>(gp.ball_centers[i]); // recall constructor
}
// assigning manifolds, with layers 100, 101, 102, 103
for (unsigned int i = 0; i < gp.n_balls; i++) {
triangulation.set_manifold ( 100 + i, balls[i] );
}
unsigned int not_curved = 1;
inner_manifold.initialize(triangulation);
triangulation.set_manifold (not_curved, inner_manifold);
triangulation.refine_global (refinement);
// ---------------- bgn setup --------------------------------------------------------------------
{
Vector<double> zeros (triangulation.n_active_cells());
std::ostringstream ss;
const std::string filename = "zeros.vtk";
DataOut<dim> data_out;
MappingQGeneric<dim> mapping (degree);
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector ( zeros, "zeros" );
data_out.build_patches (mapping, 4, DataOut<dim>::curved_inner_cells);
std::ofstream output (filename.c_str());
data_out.write_vtk (output);
}
}
template <int dim>
void Mygrid<dim>::setup_system () {
dof_handler.distribute_dofs (fe);
constraints.clear ();
DoFTools::make_zero_boundary_constraints (dof_handler, constraints);
constraints.close ();
printf ("dofs=%d,\tcells=%d\n",dof_handler.n_dofs(),triangulation.n_active_cells() );
}
// ------------------------------------ RUN ------------------------------------------------------//
template <int dim>
void Mygrid<dim>::run () {
make_grid(); //printf(" ... make_grid\n");
setup_system (); //printf(" ... setup\n");
}
int main (int argc, char **argv) {
using namespace dealii;
deallog.depth_console (0);
{
// **** Change values here (degree, refinement) ******/
Mygrid<2> problem_2d( 4,1 );
problem_2d.run ();
}
// cout << printout.str();
return 0;
}
