Hi, I have written this code for testing purpose. In case of a rotated cube, the shape functions seemed to be good as they were parallel to the edge, although there was some small round-off errors at the value 0. However, when I tested with the grid which was loaded from GMSH, the shape functions failed. In the attachment, I showed the information for the first active cell only. It is easy to get information of other cells with the codes.
I think that there exists something wrong with the cell that are not a cube, but a hexahedron when initializing the shape functions (fe_values). I think it would be easier for you to help me with the attachment below. Thank you very much. Best regards. Pham Ngoc Kien Vào Th 2, 15 thg 4, 2019 vào lúc 21:48 Wolfgang Bangerth < bange...@colostate.edu> đã viết: > On 4/11/19 9:17 PM, Phạm Ngọc Kiên wrote: > > Testing for an edge whose global vertices located from (0,0,0) to > (0,0,1) in > > real coordinates. > > With a cube I get the shape function vectors at the dof related to the > edge, > > for examples, (0,0,0), (0,0,-0.25), (0,0,-0.5), (0,0,-1), which are > parallel > > to the edge. > > However, if I create a mesh with GMSH and check the shape function > vectors, > > some of them are: (-0.0868383, -0.428, 1), (-0.129933, -0.402772, 1), > > (-0.052103 -0.2568 0.6),.., > > which are not parallel to the edge. > > > > In my limited opinion, I will get the wrong solution as the shape > functions > > created by my code are not parallel to the edge in my model. > > Could you please tell me why this happen and how to fix it? > > Good question. This does sound wrong. To figure out what is happening > exactly, > we typically construct small test cases that try to illustrate one > particular > issue. In your case, I would try to create a mesh with exactly one cell -- > either the one you have above, or maybe even simpler a cube that is just > rotated. Then create a DoFHandler on it and output the values of the shape > functions. This way, you really have only one thing that could go wrong, > and > it's easy to demonstrate what you see and how it differs from expectations. > > Do you think you could come up with a small program that does this? > > Best > WB > > -- > ------------------------------------------------------------------------ > 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.
*
* ---------------------------------------------------------------------
*
* Author: Pham Ngoc Kien, Seoul National University, 2019
*/
#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/fe/fe_nedelec.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 print_mesh_info(const Triangulation<dim> &triangulation,
const std::string &filename);
Triangulation<dim> triangulation;
MappingQ<dim> mapping;
const unsigned int order = 0; //order p of element
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
Point<dim>get_vertex_pos(const typename DoFHandler<dim>::active_cell_iterator &cell,
const FESystem<dim> &fe,
unsigned int edge,
unsigned int vertex) const;
};
template<int dim>
EM<dim>::EM()
:mapping(1),
dof_handler(triangulation),
fe(FE_NedelecSZ<dim>(order), 2)
{
}
template<int dim>
EM<dim>::~EM() {
dof_handler.clear();
}
template<int dim>
void EM<dim>::grid_generation() {
//load grid from gmsh
GridIn<dim> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("untitled.msh");
gridin.read_msh(f);
//create grid by Dealii function
//GridGenerator::hyper_cube(triangulation, 0, 1);
//GridTools::rotate(numbers::PI/4,0,triangulation);
print_mesh_info(triangulation, "model.eps");
}
template<int dim>
void EM<dim>::setup_system() {
dof_handler.distribute_dofs(fe);
}
template<int dim>
void EM<dim>::assemble_system() {
QGaussLobatto<dim> quadrature_formula(3);
FEValues<dim> fe_values(fe, quadrature_formula,
update_values | update_gradients | update_jacobians |
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;
//here I only show information in the first active cell
// if want to loop over all cells then using the for loop below
auto cell = dof_handler.begin_active();
// for (; cell != dof_handler.end(); ++cell) {
fe_values.reinit(cell);
std::cout<<"\n";
for (int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j) {
std::cout<<"cell vertex "<<j<<" has position : "<<cell->vertex(j)<<"\n";
}
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
const unsigned int block_index_i = fe.system_to_block_index(i).first;
std::cout<<"dof i_th: "<<i<<"\n";
unsigned int edge= (int)i/2.;
std::cout<<"edge number: "<<edge<<"\n";
std::cout<<"edge vertices: \n";
std::cout<<get_vertex_pos(cell,fe,edge,0)<<" and "<<get_vertex_pos(cell,fe,edge,1)<<"\n";
std::cout<<"vector from first to second vertex of the edge: "<<get_vertex_pos(cell,fe,edge,1)-get_vertex_pos(cell,fe,edge,0)<<"\n";
std::cout<<"shape functions on real cell: \n";
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
std::cout<<fe_values[vec[block_index_i]].value(i, q_point)<<"\n";
}
std::cout<<"\n";
}
// }
}
template<int dim>
void EM<dim>::run() {
grid_generation();
setup_system();
assemble_system();
}
template<int dim>
void EM<dim>::print_mesh_info(const Triangulation<dim> &triangulation, const std::string &filename) {
std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << triangulation.n_active_cells() << std::endl;
// produce a graphical representation of the mesh to an output file:
std::ofstream out(filename);
GridOut grid_out;
grid_out.write_eps(triangulation, out);
std::cout << " written to " << filename << std::endl << std::endl;
std::cout << "show information and the shape functions at first active cell only:" << std::endl;
}
template<int dim>
Point <dim> EM<dim>::get_vertex_pos(const typename DoFHandler<dim>::active_cell_iterator &cell,
const FESystem<dim> &fe,
unsigned int edge,
unsigned int vertex) const {
//local vertices of edge
unsigned int v0_loc = GeometryInfo<dim>::line_to_cell_vertices(edge, 0);
unsigned int v1_loc = GeometryInfo<dim>::line_to_cell_vertices(edge, 1);
//local positions of 2 vertices
Point<dim> v0_loc_pos = GeometryInfo<dim>::unit_cell_vertex(v0_loc);
Point<dim> v1_loc_pos = GeometryInfo<dim>::unit_cell_vertex(v1_loc);
//get the real coordinates of vertices of the edge
std::vector<Point<dim>> line_loc_vertices(2);
line_loc_vertices[0] = v0_loc_pos;
line_loc_vertices[1] = v1_loc_pos;
Quadrature<dim> q(line_loc_vertices);
FEValues<dim> fe_values_q(fe, q, update_quadrature_points);
fe_values_q.reinit(cell);
//position in the real cell
return fe_values_q.get_quadrature_points()[vertex];
}
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;
}
}
untitled.msh
Description: Mesh model
