Dear User Group,
I'm currently playing around with deal.ii to solve Maxwell's Equations and
combining Tutorial pieces
I try to solve an Eigenvalue Problem (similar to the scalar-valued step-36)
in Electrodynamics, namely find the Eigenmodes of a Spherical waveguide.
The differential equation you try to solve there is the same as in the
example in this document about Nedelec's Finite Elements (nedelec.pdf
<http://www.math.chalmers.se/~stig/underv/doktorandkurs/FEM/200607/nedelec.pdf>),
namely curl(curl(E)) = omega^2*E. Now, for 2 spacial dimensions using
Nedelec Elements with two in-plane Vector components you only calculate TE
(Transversal Electric) modes, but this works fine. If I want to calculate
the TM (Transversal Magnetic) modes, I THINK (please feel free to lecture
me) I should use the differential equation for B, namely curl(curl(B)) =
omega^2*B, but for B the normal components are required to be continuous,
so I should use Raviart_Thomas Elements instead of Nedelec and use the
Boundary condition of vanishing normal components, i.e.
project_boundary_values_div_conforming(...). This is just changing 2 lines
in the Code, but gives the error below at runtime (The Code is attached. As
Error happens quite early in the make_grid_and_dofs function and the code
is very much tutorial-like, I hope it is readable).
Additional question: is there a possibilty to use 3d-Nedelec Elements on a
2D geometry, or should I use an FE_System with a Nedelec Element and a
FE_Q? In this case, It seems I have to reformulate the curl operator, which
I'm not overly keen on...
Thanks in advance for your precious time!
Felix Schwarz
--------------------------------------------------------
An error occurred in line <77> of file
</build/deal.ii-MRHzH9/deal.ii-8.4.2/source/fe/fe_poly_tensor.cc> in
function
void dealii::{anonymous}::get_face_sign_change_rt(const cell_iterator&,
unsigned int, std::vector<double>&)
The violated condition was:
f * dofs_per_face + j < face_sign.size()
The name and call sequence of the exception was:
ExcInternalError()
Additional Information:
This exception -- which is used in many places in the library -- usually
indicates that some condition which the author of the code thought must be
satisfied at a certain point in an algorithm, is not fulfilled. An example
would be that the first part of an algorithm sorts elements of an array in
ascending order, and a second part of the algorithm later encounters an an
element that is not larger than the previous one.
There is usually not very much you can do if you encounter such an
exception since it indicates an error in deal.II, not in your own program.
Try to come up with the smallest possible program that still demonstrates
the error and contact the deal.II mailing lists with it to obtain help.
Stacktrace:
-----------
#0 /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2:
#1 /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2:
dealii::FE_PolyTensor<dealii::PolynomialsRaviartThomas<2>, 2,
2>::fill_fe_values(dealii::TriaIterator<dealii::CellAccessor<2, 2> >
const&, dealii::CellSimilarity::Similarity, dealii::Quadrature<2> const&,
dealii::Mapping<2, 2> const&, dealii::Mapping<2, 2>::InternalDataBase
const&, dealii::internal::FEValues::MappingRelatedData<2, 2> const&,
dealii::FiniteElement<2, 2>::InternalDataBase const&,
dealii::internal::FEValues::FiniteElementRelatedData<2, 2>&) const
#2 /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: dealii::FEValues<2,
2>::do_reinit()
#3 /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: void
dealii::FEValues<2, 2>::reinit<dealii::DoFHandler,
false>(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2,
2>, false> > const&)
#4 /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: void
dealii::hp::FEValues<2, 2>::reinit<dealii::DoFHandler<2, 2>,
false>(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2,
2>, false> >, unsigned int, unsigned int, unsigned int)
#5 /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.8.4.2: void
dealii::VectorTools::project_boundary_values_div_conforming<2>(dealii::DoFHandler<2,
2> const&, unsigned int, dealii::Function<2, double> const&, unsigned char,
dealii::ConstraintMatrix&, dealii::Mapping<2, 2> const&)
#6 ./Microwave_Resonator:
MR_fesc::EigenvalueProblem<2>::make_grid_and_dofs()
#7 ./Microwave_Resonator: MR_fesc::EigenvalueProblem<2>::run()
#8 ./Microwave_Resonator: main
--
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.
//Try to solve the Schrödinger Equation for the spherical Harmonic Oszi
#include <deal.II/base/logstream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.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/dofs/dof_renumbering.h>
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_raviart_thomas.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/full_matrix.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/slepc_solver.h>
#include <deal.II/grid/manifold_lib.h>
#include <set>
#include <fstream>
#include <iostream>
namespace MR_fesc{
using namespace dealii;
template<int dim>
class EigenvalueProblem{
public:
EigenvalueProblem();
void run();
private:
//unsigned int num_eigs;
void make_grid_and_dofs();
void assemble_system();
unsigned int solve();
void output_results() const;
Triangulation<dim> triang;
//for TE modes:
//FE_Nedelec<dim> fe;
//for TM modes:
FE_RaviartThomas<dim> fe;
DoFHandler<dim> dof_h;
PETScWrappers::SparseMatrix H,S;
std::vector<PETScWrappers::MPI::Vector> eigenfunctions;
std::vector<double> eigenvals;
ConstraintMatrix constraints;
};
template<int dim>
EigenvalueProblem<dim>::EigenvalueProblem():
fe(1),dof_h(triang){}
template<int dim>
void EigenvalueProblem<dim>::make_grid_and_dofs(){
GridGenerator::hyper_ball(triang); //default radius 1 around 0,0
static const SphericalManifold<dim> bound;
triang.set_all_manifold_ids_on_boundary(0);
triang.set_manifold(0,bound);
triang.refine_global(3);
GridOut gridout;
std::ofstream gridofstream("grid.eps");
gridout.write_eps(triang,gridofstream);
dof_h.distribute_dofs(fe);
//renumbering, should benchmark this at some point
//DoFRenumbering::Cuthill_McKee(dof_h);
//DoFRenumbering::block_wise(dof_h);
//BOUNDARY VALUES:
//for TE Waves, we actually compute the electric field and
//demand n x E = 0
//VectorTools::project_boundary_values_curl_conforming_l2
// (dof_h,0,ZeroFunction<dim>(dim),0,constraints); //Mapping?
//for TM Waves, we compute the magnetic field B and demand
//dot(n,B) = 0
VectorTools::project_boundary_values_div_conforming
(dof_h,0,ZeroFunction<dim>(dim),0,constraints);
constraints.close();
H.reinit(dof_h.n_dofs(),
dof_h.n_dofs(),
dof_h.max_couplings_between_dofs());
S.reinit(dof_h.n_dofs(),
dof_h.n_dofs(),
dof_h.max_couplings_between_dofs());
IndexSet eigs_ind_set = dof_h.locally_owned_dofs();
eigenfunctions.resize(13);
for(unsigned int i=0; i<eigenfunctions.size();++i){
eigenfunctions[i].reinit(eigs_ind_set,MPI_COMM_WORLD);
}
eigenvals.resize(eigenfunctions.size());
}
// template<int dim>
// std::vector<double> pot_func(const std::vector<Point<dim> > &pts,double scale){
// std::vector<double> retval(pts.size());
// for(int p=0; p<pts.size();++p){
// for(int i=0; i<dim; ++i){
// retval[p]+=scale*pts[p][i]*pts[p][i];
// }
// }
// return retval;
// }
template<int dim>
void EigenvalueProblem<dim>::assemble_system(){
QGauss<dim> quadrature_formula(3);
FEValues<dim> fe_values(fe,quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dpc = fe.dofs_per_cell;
const unsigned int nqp = quadrature_formula.size();
FullMatrix<double> cell_H(dpc,dpc);
FullMatrix<double> cell_S(dpc,dpc);
std::vector<types::global_dof_index> local_dof_indices (dpc);
//std::vector<double> potential_values(nqp);
std::vector<typename internal::CurlType<dim>::type> curl_E(dpc); //curl_type is probs not defined here
std::vector<Tensor<1,dim> > val_E(dpc); //FIXME: is that the right type?
const FEValuesExtractors::Vector Efield(0);
//weak form of Maxwell's Wave-Equation:
//(rot(E),rot(v)) = omega_h (eps*E,v) for all v
for(auto cell = dof_h.begin_active(); cell!=dof_h.end(); ++cell){
fe_values.reinit(cell);
cell_H=0;
cell_S=0;
//potential_values = pot_func(fe_values.get_quadrature_points(),500);
for(unsigned int q_point=0; q_point<nqp; ++q_point){
//just for debugging, output cell_H and cell_S
//std::cout << "Cell " << q_point << "\n";
for(unsigned int k=0; k<dpc;++k){
curl_E[k] = fe_values[Efield].curl(k,q_point);
val_E[k] = fe_values[Efield].value(k,q_point);
}
for(unsigned int i=0; i<dpc; ++i){
for(unsigned int j=0; j<dpc; ++j){
cell_H(i,j)+=(curl_E[i]*curl_E[j]*
fe_values.JxW(q_point));
cell_S(i,j)+=(val_E[i]*val_E[j]*
fe_values.JxW(q_point));
//std::cout <<
}
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_H,local_dof_indices,H);
constraints.distribute_local_to_global(cell_S,local_dof_indices,S);
}
H.compress(VectorOperation::add);
S.compress(VectorOperation::add);
double min_spurious_eigenvalue = std::numeric_limits<double>::max(),
max_spurious_eigenvalue = -std::numeric_limits<double>::max();
for (unsigned int i = 0; i < dof_h.n_dofs(); ++i)
if (constraints.is_constrained(i))
{
const double ev = H(i,i)/S(i,i);
min_spurious_eigenvalue = std::min (min_spurious_eigenvalue, ev);
max_spurious_eigenvalue = std::max (max_spurious_eigenvalue, ev);
}
std::cout << " Spurious eigenvalues are all in the interval "
<< "[" << min_spurious_eigenvalue << "," << max_spurious_eigenvalue << "]"
<< std::endl;
}
template<int dim>
unsigned int EigenvalueProblem<dim>::solve(){
SolverControl solver_control(dof_h.n_dofs(),1e-9);
SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control);
eigensolver.set_which_eigenpairs(EPS_SMALLEST_REAL);
eigensolver.set_problem_type (EPS_GHEP);
eigensolver.solve(H,S,eigenvals,eigenfunctions,eigenvals.size());
for (unsigned int i=0; i<eigenfunctions.size(); ++i){
eigenfunctions[i] /= eigenfunctions[i].linfty_norm ();
}
return solver_control.last_step ();
}
template <int dim>
void EigenvalueProblem<dim>::output_results () const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_h);
for (unsigned int i=0; i<eigenfunctions.size(); ++i)
data_out.add_data_vector (eigenfunctions[i],
std::string("eigenfunction_") +
Utilities::int_to_string(i));
data_out.build_patches ();
std::ofstream output ("eigenvectors.vtk");
data_out.write_vtk (output);
}
template<int dim>
void EigenvalueProblem<dim>::run(){
make_grid_and_dofs();
assemble_system();
const unsigned int n_iterations = solve ();
std::cout << " Solver converged in " << n_iterations
<< " iterations." << std::endl;
output_results ();
std::cout << std::endl;
for (unsigned int i=0; i<eigenvals.size(); ++i)
std::cout << " Eigenvalue " << i
<< " : " << eigenvals[i]
<< std::endl;
}
}
int main(int argc, char **argv){
using namespace MR_fesc;
using namespace dealii;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
EigenvalueProblem<2> prob;
prob.run();
return 0;
}
