Hi,
The attachment is my small program that I created.
My fe system have 24 dofs for both real and imaginary of 12 edges of a cell.
No matter how I change the parameter "order" the number of dofs is
unchanged.
I think I have problem when constructing the object.
I would like to thank you very much.
Best,
Kien
Vào Th 5, 21 thg 2, 2019 vào lúc 23:38 Wolfgang Bangerth <
bange...@colostate.edu> đã viết:
> On 2/20/19 7:01 PM, Phạm Ngọc Kiên wrote:
> > I tested in my codes using
> >
> > template<int dim>
> > CSEM<dim>::CSEM()
> > ://mapping(1),
> > dof_handler(triangulation),
> > fe(FE_NedelecSZ<dim>(order),1,//(order), multiplicity
> > FE_NedelecSZ<dim>(order),1) {}
> >
> > When I changed the parameter order, I saw only the number of dofs was
> only of lowest order dofs no matter how big the order was
> >
> > I think that there should be a way to get higher order dofs, but I
> cannot find it now.
>
> Can you create a small program that demonstrates this issue?
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
> --
> 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.
>
--
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) 2013 - 2018 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.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*/
#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/manifold_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.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/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_nedelec_sz.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/data_out_dof_data.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/fe/mapping_q.h>
#include <complex>
#include <iostream>
#include <fstream>
#include <map>
#include <complex>
#include <cmath>
using namespace dealii;
template<int dim>
class EM {
public:
EM();
~EM();
void run();
private:
void grid_generation();
void setup_system();
void assemble_system();
void solve();
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
const unsigned int order = 3; //order of edge-based element
};
template<int dim>
EM<dim>::EM()
:dof_handler(triangulation),
fe(FE_NedelecSZ<dim>(order), 1, //(order), multiplicity
FE_NedelecSZ<dim>(order), 1) {}
template<int dim>
EM<dim>::~EM() {
dof_handler.clear();
}
template<int dim>
void EM<dim>::grid_generation() {
GridGenerator::hyper_cube(triangulation, -4, 4);
}
template<int dim>
void EM<dim>::setup_system() {
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::project_boundary_values_curl_conforming_l2(dof_handler,
0,
ConstantFunction<dim>(0., fe.n_components()),
1,
constraints);
VectorTools::project_boundary_values_curl_conforming_l2(dof_handler,
dim,
ConstantFunction<dim>(0., fe.n_components()),
1,
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
}
template<int dim>
void EM<dim>::assemble_system() {
QGauss<dim> quadrature_formula(order + 2);
FEValues<dim> fe_values(fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
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);
const FEValuesExtractors::Vector Er(0);
const FEValuesExtractors::Vector Ei(dim);
std::vector<FEValuesExtractors::Vector> vec(2);
vec[0] = Er;
vec[1] = Ei;
auto cell = dof_handler.begin_active();
for (; cell != dof_handler.end(); ++cell) {
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
// block index 0 - real, 1 - imaginary
const unsigned int block_index_i = fe.system_to_block_index(i).first;
for (unsigned int j = i; j < dofs_per_cell; ++j) {
const unsigned int block_index_j = fe.system_to_block_index(j).first;
double mass_part = 0;
double curl_part = 0;
const double mu_term = 1.;
double kappa_term = 1.;
if (block_index_i == block_index_j) {
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
curl_part
+= fe_values[vec[block_index_i]].curl(i, q_point)
* fe_values[vec[block_index_j]].curl(j, q_point)
* fe_values.JxW(q_point);
}
cell_matrix(i, j) = mu_term * curl_part;
} else
{
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
mass_part
+= fe_values[vec[block_index_i]].value(i, q_point)
* fe_values[vec[block_index_j]].value(j, q_point)
* fe_values.JxW(q_point);
}
cell_matrix(i, j) = kappa_term * mass_part;
}
}
}
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
// block index 0 - real, 1 - imaginary
const unsigned int block_index_i = fe.system_to_block_index(i).first;
double rhs_term = 0;
Tensor<1, dim> rhs_value;
for (unsigned int d = 0; d < dim; ++d) {
rhs_value[d] = 0.;
}
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
rhs_term += rhs_value * fe_values[vec[block_index_i]].value(i, q_point) * fe_values.JxW(q_point);
}
cell_rhs(i) = rhs_term;
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix, cell_rhs,
local_dof_indices, system_matrix, system_rhs, false);
}
}
template<int dim>
void EM<dim>::solve() {
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix);
A_direct.vmult(solution, system_rhs);
constraints.distribute(solution);
}
template<int dim>
void EM<dim>::run() {
grid_generation();
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
setup_system();
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs() << "\n \n";
assemble_system();
solve();
}
int main() {
try {
EM<3> em_3d;
em_3d.run();
}
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;
}
}
