Dear Daniel,
Thank you for a fast answer. I think I got the matrix and right-hand side
correct now, please see test.cc for the example code.
What I want to from the start was to get the vector-valued gradient from a
scalar field, but I could not get this to work. See section "Do not work"
in function setup_system in test.cc, this gives 0 everywhere. I want to use
the information from test vector but since it was made with a different
dof_handler it doesn´t seem to work. Is there a way around this problem?
Can a vector be transferred from one dof_handler to another? or is there a
better way to solve this?
Best,
Joel
On Friday, August 12, 2016 at 12:00:36 AM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Does the equation you want to solve for have multiple components?
> Otherwise it would not make sense to multiply in the right hand side
> something with a gradient.
> If you want to project the gradient into a vector valued finite element
> ansatz space, then the matrix you are assembling looks weird.
> In what you are doing you are also considering the coupling between all
> components.
>
> You are calculating the gradient of your scalar finite element function
> correctly, but you need to find the correct component to use.
> For doing this you can use something like:
>
> const unsigned int component = fe.system_to_component_index(i).first;
> cell_rhs(i) += ((fe_values.shape_value (i, q_index) *
> fe_function_grad[q_index][component])*
> fe_values.JxW (q_index));
>
> If you want to assemble a mass matrix for a vector-valued finite element,
> you have also to restrict assembling for the matrix to the case
> where the component for the ith and jth ansatz function are the same.
>
> You might want to have a look at step-8 [1] and the module "Handling
> vector valued problems"[2].
>
> Best,
> Daniel
>
> [1]
> https://www.dealii.org/8.4.1/doxygen/deal.II/step_8.html#ElasticProblemassemble_system
> [2]
> https://www.dealii.org/8.4.1/doxygen/deal.II/group__vector__valued.html
>
--
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) 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.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 1999
*/
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_accessor.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/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_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/compressed_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <deal.II/base/logstream.h>
using namespace dealii;
const double pi=3.14159265359;
double center_id=0;
template <int dim>
class Step4
{
public:
Step4 ();
void run ();
private:
void make_grid ();
void setup_system();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation<dim> triangulation;
FE_Q<dim> fe_scalar;
FESystem<dim> fe;
DoFHandler<dim> dof_handler_scalar;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
ConstraintMatrix constraints;
Vector<double> solution;
Vector<double> system_rhs;
Vector<double> test;
std::vector<Point<dim>> nodes;
};
template <int dim>
Step4<dim>::Step4 ()
:
fe_scalar (1),
fe (FE_Q<dim>(1),dim),
dof_handler (triangulation),
dof_handler_scalar (triangulation)
{}
template <int dim>
void Step4<dim>::make_grid ()
{
GridGenerator::hyper_cube (triangulation, -2.5, 2.5);
triangulation.refine_global (4);
for (unsigned int step=0; step<2; ++step)
{
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int v=0;v < GeometryInfo<dim>::vertices_per_cell;++v)
{
Point<dim> center;
for (unsigned int i=0; i<dim; ++i)
center(i)=0;
const double distance_from_center = center.distance (cell->vertex(v));
if (std::fabs(distance_from_center) < 1.0)
{
cell->set_refine_flag ();
break;
}
}
}
triangulation.execute_coarsening_and_refinement ();
}
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
}
template <int dim>
void Step4<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
dof_handler_scalar.distribute_dofs (fe_scalar);
std::cout << " Number of degrees of freedom for vector: "
<< dof_handler.n_dofs()
<< std::endl;
std::cout << " Number of degrees of freedom for scalar: "
<< dof_handler_scalar.n_dofs()
<< std::endl;
constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,constraints);
constraints.close ();
CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, c_sparsity,constraints,/*keep_constrained_dofs = */false);
sparsity_pattern.copy_from(c_sparsity);
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
/*
// Do not work
test.reinit (dof_handler_scalar.n_dofs());
nodes.resize(dof_handler_scalar.n_dofs());
MappingQ1<dim> mapping;
DoFTools::map_dofs_to_support_points(mapping,dof_handler_scalar,nodes);
*/
// Works
test.reinit (dof_handler.n_dofs());
nodes.resize(dof_handler.n_dofs());
MappingQ1<dim> mapping;
DoFTools::map_dofs_to_support_points(mapping,dof_handler,nodes);
for (unsigned int i=0;i<test.size();i++)
{
Point<dim> center;
for (unsigned int i=0; i<dim; ++i)
center(i)=0;
if (nodes[i].distance(center) < 1e-8)
{
center_id = i;
std::cout << "center id: " << center_id << std::endl;
break;
}
}
std::cout << "fe_degree: " << dof_handler.get_fe().degree << std::endl;
std::cout << "quad_degree: " << dof_handler.get_fe().degree+1 << std::endl;
}
template <int dim>
void Step4<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(dof_handler.get_fe().degree+1);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
std::vector<Tensor<1, dim>> fe_function_grad(quadrature_formula.size());
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);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
fe_values.get_function_gradients(test,fe_function_grad);
cell_matrix = 0;
cell_rhs = 0;
const FEValuesExtractors::Vector velocities (0);
for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q_index);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> phi_j_u = fe_values[velocities].value (j, q_index);
cell_matrix(i,j) += ((phi_i_u*phi_j_u)*
fe_values.JxW (q_index));
}
cell_rhs(i) += ((phi_i_u*fe_function_grad[q_index])*
fe_values.JxW (q_index));
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
template <int dim>
void Step4<dim>::solve ()
{
SolverControl solver_control (1000, 1e-6);
SolverCG<> solver (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
solver.solve (system_matrix, solution, system_rhs, preconditioner);
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;
constraints.distribute (solution);
}
template <int dim>
void Step4<dim>::output_results () const
{
// solution
DataOut<dim> data_out;
// vector
std::vector<std::string> solution_names;
solution_names.push_back ("x");
solution_names.push_back ("y");
solution_names.push_back ("z");
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names);
data_out.build_patches ();
std::ofstream output ("./solution.vtk");
data_out.write_vtk (output);
//test
DataOut<dim> data_out2;
data_out2.attach_dof_handler (dof_handler);
data_out2.add_data_vector (test, "solution");
data_out2.build_patches ();
std::ofstream output2 ("./test.vtk");
data_out2.write_vtk (output2);
//text file
std::ofstream myfile;
myfile.open ("./output/output.txt", std::ios::out | std::ios::trunc);
for (unsigned int i=0;i<dof_handler.n_dofs();i++)
{
if (fabs(nodes[i](1)) < 1e-8 and fabs(nodes[i](2)) < 1e-8)
{
myfile << nodes[i] << "\t" << solution[i] << std::endl;
}
}
myfile.close();
}
template <int dim>
void Step4<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
make_grid();
setup_system ();
for (unsigned int i=0;i<test.size();i++)
{
Point<dim> center;
for (unsigned int i=0; i<dim; ++i)
center(i)=0;
//Gaussian test function
test(i)= exp(-pow((nodes[i].distance(center)),2.0));
}
assemble_system ();
solve ();
output_results ();
}
int main ()
{
deallog.depth_console (0);
Step4<3> test;
test.run();
return 0;
}
