Some context:
I am trying to model a waveguide which would require me to solve the double
curl equation curlxcurlxE-cE=0 where c is some constant(real or imaginary)
subject to nxE=nxf on different parts of the boundary where f can be real
or imaginary also.
The problem requires the use of Nedelec elements but you already know that.
Since the problem is complex I am breaking it down into two with FEsystem.
I am trying to approach this by solving this simple system before
complicating things:
curlxcurlxE_r-cE_r=0
nxE_r=1 on an "inlet"
nxE_r=0 everywhere else
and
curlxcurlxE_i-cE_i=0
nxE_i=1 on an "inlet"
nxE_i=0 everywhere else
These two systems are uncoupled and exactly the same. The issue has been
that my BCs aren't being applied for the second part of the equation. I am
using VectorTools::project_boundary_values_curl_conforming_l2(). The first
part of the equation gives a solution (inaccurate but that's not my issue
here) with the BCs applied.
I print out the sparsity pattern and there seems to be a difference between
applying BCs and not applying them which makes me conclude that
project_boundary_values_curl_conforming_l2 is effecting the matrix.
I printed out the global right hand side for a smaller problem and the
second half of it, which is the half that is related to the second problem,
is just zeros so my guess is that
constraints.distribute_local_to_global(cell_matrix, cell_rhs,
local_dof_indices, system_matrix, system_rhs) isn't applying the BC in the
cell_rhs when FENedelec systems are involved.
I have attached my code here.
Thanks
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/8c95e3cb-1643-4cbf-a79f-b7ebc7f02104n%40googlegroups.com.
#include <deal.II/base/function.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h>
#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_out.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <iostream>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <vector>
using namespace dealii;
// Perfect conductor BC nxE=0
class PerfConductor : public Function<3>
{
public:
PerfConductor(): Function<3>(3)
{}
void vector_value_list(const std::vector<Point<3>> &points,
std::vector<Vector<double>> & values) const override;
};
void PerfConductor::vector_value_list(const std::vector<Point<3>> &points,
std::vector<Vector<double>> & values) const
{
for (unsigned int i = 0; i < points.size(); ++i)
{
values[i][0] = 0.0;
values[i][1] = 0.0;
values[i][2] = 0.0;
}
}
//Inlet BC nXE=1
class Inlet: public Function<3>
{
public:
Inlet(): Function<3>(3)
{}
void vector_value_list(const std::vector<Point<3>> &points,
std::vector<Vector<double>> & values) const override;
};
void Inlet::vector_value_list(const std::vector<Point<3>> &points,
std::vector<Vector<double>> & values) const
{
for (unsigned int i = 0; i < points.size(); ++i)
{
values[i][0] = 1.0;
values[i][1] = 0.0;
values[i][2] = 0.0;
}
}
class Nedlec
{
public:
Nedlec();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<3> triangulation;
DoFHandler<3> dof_handler;
FESystem<3> fe;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
//Constructor
Nedlec::Nedlec()
: dof_handler(triangulation)
, fe(FE_Nedelec<3>(0), 2)
{}
//makegrid
void Nedlec::make_grid()
{
Point<3> p1(0,0,0) , p2 (0.2,0.1,1) ;
GridGenerator::hyper_rectangle(triangulation,p1,p2);
triangulation.begin_active()->face(4)->set_boundary_id(1); //set the 4th boundary id to 1 to apply dirichlet BC
triangulation.refine_global(2);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell->set_refine_flag(RefinementCase<3>::cut_z);
}
triangulation.execute_coarsening_and_refinement();
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell->set_refine_flag(RefinementCase<3>::cut_z);
}
triangulation.execute_coarsening_and_refinement();
std::ofstream out("grid-1.vtk");
GridOut grid_out;
grid_out.write_vtk(triangulation, out);
std::cout << "Grid written to grid-1.vtk" << std::endl;
std::cout << "Number of active cells = " << triangulation.n_active_cells() << std::endl ;
}
void Nedlec::setup_system()
{
dof_handler.distribute_dofs(fe);
constraints.clear();
PerfConductor perfBC;
//Apply perfect conductor BC on the walls with ID0 and for the real part of the eqaution
VectorTools::project_boundary_values_curl_conforming_l2(
dof_handler,
0,
perfBC,
0,
constraints,
StaticMappingQ1<3>::mapping);
//Apply perfect conductor BC on the walls with ID0 and for the imaginary part of the eqaution
VectorTools::project_boundary_values_curl_conforming_l2(
dof_handler,
3,
perfBC,
0,
constraints,
StaticMappingQ1<3>::mapping);
std::cout << "PEC was applied" << std::endl;
Inlet InletBC;
//Apply inlet BC on the walls with ID1 and for the real part of the eqaution
VectorTools::project_boundary_values_curl_conforming_l2(
dof_handler,
0,
InletBC,
1,
constraints,
StaticMappingQ1<3>::mapping);
//Apply inlet BC on the walls with ID1 and for the complex part of the eqaution
VectorTools::project_boundary_values_curl_conforming_l2(
dof_handler,
3,
InletBC,
1,
constraints,
StaticMappingQ1<3>::mapping);
std::cout << "Inlet BC was applied" << std::endl;
constraints.close();
DynamicSparsityPattern sp(dof_handler.n_dofs(),dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
sp,
constraints,
false);
sparsity_pattern.copy_from(sp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
//can only print psarsity pattern for small meshes. St golbal refinement to 2 when using this
std::ofstream out("sparsity_pattern.svg");
sparsity_pattern.print_svg(out);
std::cout << "Number of degress of freedom = " << dof_handler.n_dofs() << std::endl ;
}
void Nedlec::assemble_system()
{
QGauss<3> quad(1*1+1); //Recommended qguass but not sure
FEValues<3> fe_values(fe, quad,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values);
//extractors for real and imaginary part
const FEValuesExtractors::Vector realpart(0);
const FEValuesExtractors::Vector imagpart(3);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
std::cout<< "Dofs per cell = " << dofs_per_cell << std::endl;
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
for (const unsigned int q: fe_values.quadrature_point_indices())
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
cell_matrix(i, j) +=
(fe_values[realpart].curl(i, q) * fe_values[realpart].curl(j, q) //curl u curl v
-
1757*fe_values[realpart].value(i, q) * fe_values[realpart].value(j, q) //u v
+
fe_values[imagpart].curl(i, q) * fe_values[imagpart].curl(j, q) //curl u curl v
-
1757*fe_values[imagpart].value(i, q) * fe_values[imagpart].value(j, q) //u v
)
*
fe_values.JxW(q) ;
}
}
}
cell->get_dof_indices(local_dof_indices);
//Apply BC
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
std::cout<<"System Assemble complete" << std::endl;
}
void Nedlec::solve()
{ std::cout<< "Solving" << std::endl ;
SparseDirectUMFPACK direct;
direct.initialize(system_matrix);
direct.vmult(solution, system_rhs);
constraints.distribute(solution);
std::cout << "Soution ended" <<std::endl;
}
void Nedlec::output_results() const
{
std::vector<std::string> solution_names(3, "E_real");
for(int i=0; i<3; i++)
solution_names.emplace_back("E_imag");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
interpretation(6,
DataComponentInterpretation::component_is_part_of_vector);
interpretation[3]=(DataComponentInterpretation::component_is_part_of_vector);
DataOut<3> data_out;
data_out.add_data_vector(dof_handler,
solution,
solution_names,
interpretation);
data_out.build_patches();
std::ofstream output("solution.vtu");
data_out.write_vtu(output);
// DataOut<3> data_out;
// data_out.attach_dof_handler(dof_handler);
// data_out.add_data_vector(solution, "solution");
// data_out.build_patches();
// std::ofstream output("solution.vtk");
// data_out.write_vtk(output);
}
void Nedlec::run()
{
make_grid() ;
setup_system() ;
assemble_system();
solve();
output_results();
}
int main()
{
Nedlec my_nedlec;
my_nedlec.run();
}
