Dear Jean-Paul,
Thank you for your reply,
I already set up the Identity preconditioner with a block matrix, but the
problem remains the same. I am not sure this is the problem because I
didn't need to set up the identity operator to any matrix in the scenario
where I could invert. I think my problem is in the Schur preconditioner
operator before to invert.
I attach the problematic code and the "working code" I am coding. The only
difference between them is the line 1251:
(const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S
); //(problematic case)
(const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S
); //(Case I could invert)
(Please excuse me for the extension of the code, but the solver function
where I have the issue is too short. In the case you consider it necessary,
I could code a shorter one)
Thank you so much,
Regards,
Felipe
El miércoles, 16 de junio de 2021 a las 15:45:27 UTC+8, Jean-Paul Pelteret
escribió:
> Dear Felipe,
>
> It might be that you need to set up the preconditioner operator with an
> exemplar matrix (since the identity operator doesn't know the size of the
> range and domain that it's working on).
>
> If that's not the issue then could you please try to reproduce this as a
> minimal example and share it with us? The program doesn't need to produce
> any meaningful result, but it would be good if it shows both the working
> scenario and the problematic one. That way we could investigate further and
> try to figure out what's going wrong here.
>
> Best,
> Jean-Paul
>
> Sent from my mobile device. Please excuse my brevity and any typos.
>
> On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo, <[email protected]>
> wrote:
>
>> Dear Jean-Paul,
>>
>> Thank you for your reply and the complete information you provide me.
>> Effectively, declaring the preconditioner (out-of-line) was one problem,
>> but the inverse_operator function is still not matching.
>>
>> What I just realized is that when I compute the operator as a
>> multiplication of a linear operator, an inverse operator and a linear
>> operator ( for instance, when I compute the Schur complement )
>>
>> const auto op_S = op_BT * op_M_inv * op_B;
>>
>> this operator becomes to the type
>> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
>>
>> ,
>>
>> but when I compute the operator as a multiplication of other linear
>> operators directly (for instance, when I compute the Schur preconditioner,
>> following the procedure in step 20 )
>>
>> const auto op_pre = linear_operator<LA::MPI::Vector>(prec_M);
>> const auto op_aS = op_BT * op_pre * op_B;
>>
>> this operator becomes to the type
>> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
>>
>> So, If I use the inverse_operator function with the first operator
>> (op_S), I can obtain the inverse without any problem,
>> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S
>> );
>>
>> But if I do it with the second case, it doesn't work because the
>> functions are not matching.
>> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S
>> );
>>
>> I am not sure what is happening because all the operators are declared
>> as "LA::MPI::Vector".
>> If you have any suggestion, would be greatly appreciated.
>>
>> Thank you so much,
>>
>> Regards,
>> Felipe Giraldo
>>
>>
>>
>> El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret
>> escribió:
>>
>>> Hi again Feilpe,
>>>
>>> Regarding the lack of documentation, I’ve opened an issue on Github to
>>> track this. You can find that here:
>>> https://github.com/dealii/dealii/issues/12466
>>>
>>> Best,
>>> Jean-Paul
>>>
>>> On 15. Jun 2021, at 12:57, Jean-Paul Pelteret <[email protected]>
>>> wrote:
>>>
>>> Hi Feilpe,
>>>
>>> Firstly, I agree that the documentation is very light on details on how
>>> to use the linear operators with Trilinos linear algebra types. We can
>>> definitely improve this, and any contributions to that effect would be
>>> greatly appreciated!
>>>
>>> Right now, I can direct you to a few of the tests that use Trilinos LA
>>> in conjunction with the inverse operator, so that you can compare what you
>>> have and figure out what the problematic differences are. There is this
>>> one, for instance
>>>
>>> https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
>>> that looks like a similar setup to yours, and
>>>
>>> https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
>>> that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).
>>>
>>> Comparing to both of these, I think that the important point might be
>>> that the preconditioner must be declared (out-of-line) before the
>>> inverse_operation() function is called. The linear operators typically
>>> expect the lifetime of the LA objects to exceed that of the operators, and
>>> so you have to first create the matrix or preconditioner and then pass it
>>> to the linear operators. This is similar to what you’d done when setting up
>>> op_M,
>>> for example. The case of passing in a deal::PreconditionerIdentity()
>>> for serial operations is a special case, and I don’t think that we’ve
>>> duplicated that for TrilinosWrappers:: PreconditionerIdentity(). Maybe
>>> that could be improved too.
>>>
>>> I hope that with this single change you’d be able to get your
>>> program running. If not, then please do let us know so that we can try to
>>> help further.
>>>
>>> Best,
>>> Jean-Paul
>>>
>>>
>>> On 14. Jun 2021, at 19:18, Juan Felipe Giraldo <[email protected]>
>>> wrote:
>>>
>>> Hello everyone,
>>>
>>> I have implemented a residual-minimization framework that somehow is
>>> similar to DPG. I want to extend my results by using parallelization using
>>> MPI with PETSc or Trilinos.
>>> So far, I have solved the saddle point problem using the Schur
>>> complement exactly how it is described in step 20. Now, I am trying to
>>> replicate exactly the same solver but using the MPI wrappers and the linear
>>> operators.
>>>
>>> The problem is that when I am trying to implement the inverse_operator
>>> to compute the Preconditioner of the Schur complement, I get an error
>>> saying that the functions are not matching "inverse_operator(op_aS,
>>> solver_aS, TrilinosWrappers::PreconditionIdentity())."
>>>
>>> There is no much documentation about linear operators in parallel
>>> solvers, so if anyone has any suggestion on how to fix this problem, it
>>> would be well appreciated.
>>>
>>> I have pasted the complete function in below:
>>>
>>>
>>> template <int dim>
>>> void FEMwDG<dim>::
>>> solve()
>>> {
>>> TimerOutput::Scope t(computing_timer, "solve");
>>>
>>> LA::MPI::PreconditionJacobi prec_M;
>>> LA::MPI::PreconditionJacobi::AdditionalData data;
>>> prec_M.initialize(system_matrix.block(0, 0), data);
>>> }
>>>
>>> auto &E = solution.block(0);
>>> auto &U = solution.block(1);
>>> const auto &L = system_rhs.block(0);
>>> const auto M =
>>> linear_operator<LA::MPI::Vector>(system_matrix.block(0, 0));
>>>
>>> const auto op_M = linear_operator(M, prec_M);
>>> const auto op_B =
>>> linear_operator<LA::MPI::Vector>(system_matrix.block(0, 1));
>>> const auto op_BT =
>>> linear_operator<LA::MPI::Vector>(system_matrix.block(1, 0));
>>>
>>> ReductionControl inner_solver_control2(5000,
>>> 1e-18 *
>>> system_rhs.l2_norm(),
>>> 1.e-2);
>>>
>>> SolverCG<LA::MPI::Vector> cg(inner_solver_control2);
>>> const auto op_M_inv = inverse_operator(M, cg, prec_M);
>>>
>>> const auto op_S = op_BT * op_M_inv * op_B;
>>> const auto op_aS = op_BT * linear_operator<LA::MPI::Vector>(prec_M)
>>> * op_B;
>>>
>>> IterationNumberControl iteration_number_control_aS(30, 1.e-18);
>>> SolverCG<LA::MPI::Vector> solver_aS(iteration_number_control_aS);
>>>
>>> const auto preconditioner_S =
>>> inverse_operator(op_aS, solver_aS,
>>> TrilinosWrappers::PreconditionIdentity());
>>> const auto schur_rhs = op_BT * op_M_inv * L ;
>>>
>>> SolverControl solver_control_S(2000, 1.e-12);
>>> SolverCG<LA::MPI::Vector> solver_S(solver_control_S);
>>>
>>> const auto op_S_inv = inverse_operator(op_S, solver_S,
>>> preconditioner_S );
>>>
>>> U = op_S_inv * schur_rhs;
>>> std::cout << solver_control_S.last_step()
>>> << " CG Schur complement iterations to obtain convergence."
>>> << std::endl;
>>> E = op_M_inv * (L - op_B * U);
>>> }
>>>
>>> Thank you so much,
>>>
>>> Regards,
>>> Felipe
>>>
>>>
>>>
>>>
>>>
>>> --
>>> 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/40b8db07-5243-4603-825b-9e7850580206n%40googlegroups.com
>>>
>>> <https://groups.google.com/d/msgid/dealii/40b8db07-5243-4603-825b-9e7850580206n%40googlegroups.com?utm_medium=email&utm_source=footer>
>>> .
>>>
>>>
>>>
>>> --
>> 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/2ece64b7-cad5-4b2b-bf61-442c7cbb2ff4n%40googlegroups.com
>>
>> <https://groups.google.com/d/msgid/dealii/2ece64b7-cad5-4b2b-bf61-442c7cbb2ff4n%40googlegroups.com?utm_medium=email&utm_source=footer>
>> .
>>
>
--
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/e1f9ce07-0c22-4af4-b302-d640ce5cad14n%40googlegroups.com.
//ASFEM method by Juan Felipe Giraldo
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/generic_linear_algebra.h>
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
#include <deal.II/lac/block_linear_operator.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/linear_operator_tools.h>
#include <deal.II/lac/packaged_operation.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include<deal.II/lac/schur_complement.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.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_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <fstream>
#include <iostream>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.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/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/assembler.h>
#include <deal.II/meshworker/loop.h>
//#include "Functions2.cc"
#include <deal.II/base/function.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/lac/vector.h>
#include <cmath>
using namespace std;
using namespace dealii;
// ---------------------Beta ------------------------------
template <int dim>
Tensor<1,dim> beta (const Point<dim>)
{
Point<dim> field;
field(0) = 1;
field(1) = 0;
return field;
}
// ---------------------Kappa ------------------------------
template <int dim>
double coefficient(const Point<dim> &p)
{
return 1E-1;
}
template <int dim>
class RHS : public Function<dim>
{
public:
RHS() : Function<dim>(1)
{}
virtual double value(const Point<dim> &p,
const unsigned int component = 0 ) const;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues() : Function<dim>(1)
{}
virtual double value(const Point<dim> &p,
const unsigned int component = 0 ) const;
};
//----------------------------RHS -----------------------------------
template <int dim>
double
RHS<dim>::
value(const Point<dim> &p,
const unsigned int ) const
{
return 0.;
}
//-------- ----------Drichlet Boundary-------------------------------
template <int dim>
double BoundaryValues<dim>:: value(const Point<dim> &p,
const unsigned int ) const
{
const double x = p[0];
const double y = p[1];
//Eriksson
if (x < 1E-15)
{
const double kappa = coefficient<dim>(p);
double r1 = (1. + sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double r2 = (1. - sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
return (exp(r1*(x-1.)) - exp(r2*(x-1.)))/(exp(-r1)-exp(-r2)) * sin(M_PI*y);
}
else
{
return 0;
}
}
// -------------Exact Solution-------------------------------
template <int dim>
class ExactSolution : public Function<dim>
{
public:
ExactSolution() : Function<dim>() //first for the error, second for the scalar u value
{}
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override;
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
void ExactSolution<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
const double x = p[0];
const double y = p[1];
values(0) = 0; //error "exact solution"
const double kappa = coefficient<dim>(p);
double r1 = (1. + sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double r2 = (1. - sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double val = (exp(r1*(x-1.)) - exp(r2*(x-1.)))/(exp(-r1)-exp(-r2)) * sin(M_PI*y);
values(1) = val;
}
template <int dim>
double ExactSolution<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
const double x = p[0];
const double y = p[1];
const double kappa = coefficient<dim>(p);
double r1 = (1. + sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double r2 = (1. - sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double val = (exp(r1*(x-1.)) - exp(r2*(x-1.)))/(exp(-r1)-exp(-r2)) * sin(M_PI*y);
return val;
}
//-----------------------------------CLASS TEMPLATE--FEMwDG-------------------------------------------
template <int dim>
class FEMwDG
{
public:
FEMwDG (const unsigned int degree);
void run (const unsigned int refinement_initial, const unsigned int typecase , const unsigned int n_cycles, const unsigned int refine_type);
~FEMwDG();
void run();
private:
using VectorType = TrilinosWrappers::MPI::Vector;
using MatrixType = TrilinosWrappers::SparseMatrix;
void make_grid(const unsigned int refinment_initial ,
const unsigned int cycle ,
const unsigned int refine_type);
void setup_system();
void exact_plot(const unsigned int cycle);
void assemble_system(const unsigned int typecase);
void assemble_cell_terms(const FEValues<dim> &cell_fe,
FullMatrix<double> &cell_matrix,
Vector<double> &cell_vector,
const unsigned int &typecase,
const double & h);
void assemble_boundary_terms(const FEFaceValues<dim> &face_fe,
FullMatrix<double> &local_matrix,
Vector<double> &local_vector,
const double &h);
//const unsigned int &typecase,
void assemble_face_terms(const FEFaceValuesBase<dim> &fe_face_values,
const FEFaceValuesBase<dim> &fe_neighbor_face_values,
FullMatrix<double> &vi_ui_matrix,
FullMatrix<double> &vi_ue_matrix,
FullMatrix<double> &ve_ui_matrix,
FullMatrix<double> &ve_ue_matrix,
const unsigned int &typecase,
const double & h);
void estimate(const unsigned int typecase);
void distribute_local_flux_to_global(
const FullMatrix<double> & vi_ui_matrix,
const FullMatrix<double> & vi_ue_matrix,
const FullMatrix<double> & ve_ui_matrix,
const FullMatrix<double> & ve_ue_matrix,
const std::vector<types::global_dof_index> & local_dof_indices,
const std::vector<types::global_dof_index> & local_neighbor_dof_indices);
void solve();
//void estimate(const unsigned int typecase);
void output_results(const unsigned int cycle) const;
const unsigned int degree;
//const unsignet int typecase;
parallel::distributed::Triangulation<dim> triangulation;
FESystem<dim> fe;
FE_DGQ<dim> fe_e;
FE_Q<dim> fe_u;
DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_e;
DoFHandler<dim> dof_handler_u;
LA::MPI::BlockSparseMatrix system_matrix;
LA::MPI::BlockSparseMatrix preconditioner_matrix;
LA::MPI::BlockVector locally_relevant_solution;
LA::MPI::BlockVector system_rhs;
LA::MPI::BlockVector solution;
std::vector<IndexSet> owned_partitioning;
std::vector<IndexSet> relevant_partitioning;
AffineConstraints<double> constraints;
AffineConstraints<double> constraints_u;
Vector<double> estimated_error_per_cell;
//BlockVector<double> estimates;
Vector<double> phi_0;
ConditionalOStream pcout;
TimerOutput computing_timer;
SolverControl solver_control;
TrilinosWrappers::SolverDirect solver;
const RHS<dim> rhs_function;
const BoundaryValues<dim> boundary_function;
const ExactSolution<dim> exact_solution;
};
template <int dim>
FEMwDG<dim>::
FEMwDG(const unsigned int degree)
:
degree(degree),
//n_refine(n_refine),
triangulation(MPI_COMM_WORLD,
typename Triangulation<dim>::MeshSmoothing
(Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening)),
fe(FE_DGQ<dim>(degree), 1,FE_Q<dim>(degree), 1), // test(e)->(D.G)- trial(u)->(C.G). Second term refers to a number of copies of fe, copy for direction (in this case is 1).
//fe(FE_DGQ<dim>(degree), 1,FE_DGQ<dim>(degree), 1), // --> ONLY DG
fe_e(degree),
fe_u(degree),
//fe(FESystem<dim>(FE_DGQ<dim>(degree), dim) , 1,FE_DGQ<dim>(degree),1),
dof_handler(triangulation),
dof_handler_e (triangulation),
dof_handler_u (triangulation),
pcout(std::cout,
Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
computing_timer(MPI_COMM_WORLD,
pcout,
TimerOutput::summary,
TimerOutput::wall_times),
solver_control(1),
solver(solver_control),
rhs_function(),
boundary_function()
//estimates(1)
{}
template <int dim>
FEMwDG<dim>::
~FEMwDG()
{
dof_handler.clear();
}
//-----------------------------------------MESH CREATION --------------------------
template <int dim>
void FEMwDG<dim>::make_grid(const unsigned int refinement_initial , const unsigned int cycle, const unsigned int refine_type)
{
if (cycle==0)
{
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(refinement_initial);
}
else
{
if (refine_type==0) // Global refinement
triangulation.refine_global(1);
if (refine_type==1) // Adaptive refinement (Dorfler marking)
{triangulation.refine_global(1);
// parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
//triangulation, estimated_error_per_cell, 0.3, 0.1);
//GridRefinement::refine_and_coarsen_fixed_fraction(
//triangulation, estimates.block(0), 0.5, 0.3);
//triangulation.execute_coarsening_and_refinement();
}
}
}
//-----------------------------------------RUN CLASS--------------------------
template <int dim>
void
FEMwDG<dim>::
setup_system()
{
TimerOutput::Scope t(computing_timer, "setup");
dof_handler.distribute_dofs(fe);
DoFRenumbering::component_wise(dof_handler);
std::vector<types::global_dof_index> dofs_per_block(2); //Grados de libertad por componentes. (escalar - escalar)
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block);//, block_component);
const unsigned int n_e = dofs_per_block[0],
n_u = dofs_per_block[1];
//std::cout << "Number of active cells: " << triangulation.n_active_cells()
// << std::endl
std::cout << "Total number of cells: " << triangulation.n_active_cells()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << n_e << '+' << n_u << ')' << std::endl;
// set the configuration for the projection - exact solution
dof_handler_u.distribute_dofs(fe_u);
DynamicSparsityPattern dsp_u (n_u, n_u);
phi_0.reinit (n_u);
constraints_u.clear ();
DoFTools::make_flux_sparsity_pattern (dof_handler_u, dsp_u, constraints_u, false);
DoFRenumbering::component_wise (dof_handler_u);
owned_partitioning.resize(2);
owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, n_e);
owned_partitioning[1] = dof_handler.locally_owned_dofs().get_view(n_e, n_e + n_u);
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
relevant_partitioning.resize(2);
relevant_partitioning[0] = locally_relevant_dofs.get_view(0, n_e);
relevant_partitioning[1] = locally_relevant_dofs.get_view(n_e, n_e + n_u);
constraints.clear();
constraints.close();
system_matrix.clear();
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
DoFTools::make_flux_sparsity_pattern(dof_handler,
dsp);
SparsityTools::distribute_sparsity_pattern(
dsp,
dof_handler.locally_owned_dofs(),
MPI_COMM_WORLD,
locally_relevant_dofs);
system_matrix.reinit(owned_partitioning,
dsp,
MPI_COMM_WORLD);
preconditioner_matrix.clear();
DoFTools::make_flux_sparsity_pattern(dof_handler,
dsp);
SparsityTools::distribute_sparsity_pattern(
dsp,
Utilities::MPI::all_gather(MPI_COMM_WORLD,dof_handler.locally_owned_dofs()),
MPI_COMM_WORLD,
locally_relevant_dofs);
preconditioner_matrix.reinit(owned_partitioning,
dsp,
MPI_COMM_WORLD);
locally_relevant_solution.reinit(owned_partitioning,
relevant_partitioning,
MPI_COMM_WORLD);
system_rhs.reinit(owned_partitioning, MPI_COMM_WORLD);
solution.reinit(owned_partitioning, MPI_COMM_WORLD);
}
double compute_penalty(const unsigned int fe_degree,
const double cell_extend_left,
const double cell_extend_right)
{
const double degree = std::max(1., static_cast<double>(fe_degree));
return degree * (degree + 1.) * 0.5 *
(1. / cell_extend_left + 1. / cell_extend_right);
}
template <int dim>
void
FEMwDG<dim>::
assemble_system(const unsigned int typecase)
{
TimerOutput::Scope t(computing_timer, "assembly");
QGauss<dim> quadrature_formula(fe.degree+1);
QGauss<dim-1> face_quadrature_formula(fe.degree+1);
const UpdateFlags update_flags = update_values
| update_gradients
| update_quadrature_points
| update_JxW_values;
const UpdateFlags face_update_flags = update_values
| update_gradients
| update_normal_vectors
| update_quadrature_points
| update_JxW_values;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<types::global_dof_index> local_neighbor_dof_indices(dofs_per_cell);
FEValues<dim> fe_values(fe, quadrature_formula, update_flags);
FEFaceValues<dim> fe_face_values(fe,face_quadrature_formula,
face_update_flags);
FEFaceValues<dim> fe_neighbor_face_values(fe, face_quadrature_formula,
face_update_flags);
FESubfaceValues<dim> fe_subface_values(fe, face_quadrature_formula,
face_update_flags);
FullMatrix<double> local_matrix(dofs_per_cell,dofs_per_cell);
Vector<double> local_vector(dofs_per_cell);
FullMatrix<double> vi_ui_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> vi_ue_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> ve_ui_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> ve_ue_matrix(dofs_per_cell, dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for(; cell!=endc; cell++)
{
if(cell->is_locally_owned())
{
local_matrix = 0;
local_vector = 0;
fe_values.reinit(cell);
double h = cell->diameter();
//------------DOMAIN INTEGRATOR----------------------------------------------
assemble_cell_terms(fe_values,
local_matrix,
local_vector,
typecase,
h);
cell->get_dof_indices(local_dof_indices);
for(unsigned int face_no=0;
face_no< GeometryInfo<dim>::faces_per_cell;
face_no++)
{
//-------------EXTERNAL BOUNDARY INTEGRATOR-------------------------------------
typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
if(face->at_boundary() )
{
const double extent1 = cell->measure() / cell->face(face_no)->measure();
const double penalty = compute_penalty(degree, extent1, extent1);
fe_face_values.reinit(cell, face_no);
assemble_boundary_terms(fe_face_values,
local_matrix,
local_vector,
penalty);
//typecase,h);
}
//-------------SKELETON INTEGRATOR ----------------------------------------------
else
{
Assert(cell->neighbor(face_no).state() == IteratorState::valid,
ExcInternalError());
typename DoFHandler<dim>::cell_iterator neighbor =
cell->neighbor(face_no);
if(face->has_children())
{
const unsigned int neighbor_face_no =
cell->neighbor_of_neighbor(face_no);
for(unsigned int subface_no=0;
subface_no < face->number_of_children();
++subface_no)
{
typename DoFHandler<dim>::cell_iterator neighbor_child =
cell->neighbor_child_on_subface(face_no, subface_no);
Assert(!neighbor_child->has_children(), ExcInternalError());
vi_ui_matrix = 0;
vi_ue_matrix = 0;
ve_ui_matrix = 0;
ve_ue_matrix = 0;
fe_subface_values.reinit(cell, face_no, subface_no);
fe_neighbor_face_values.reinit(neighbor_child,
neighbor_face_no);
const double extent1 = cell-> measure() / cell-> face(face_no)-> measure();
const double extent2 = neighbor_child -> measure() / neighbor_child -> face(subface_no)-> measure();
const double penalty = compute_penalty(degree, extent1, extent2);
assemble_face_terms(fe_subface_values,
fe_neighbor_face_values,
vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
typecase,
penalty);
//h);
neighbor_child->get_dof_indices(local_neighbor_dof_indices);
distribute_local_flux_to_global(vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
local_dof_indices,
local_neighbor_dof_indices);
}
}
else
{
//---------------ENSAMBLANDO CONTRIBUCIÓN LOCAL EN LOS BORDES AL SISTEMA GLOBAL CUANDO EXISTE VECINO-------------------------------------
if(neighbor->level() == cell->level() &&
neighbor->index() > cell->index() )
{
const unsigned int neighbor_face_no =
cell->neighbor_of_neighbor(face_no);
vi_ui_matrix = 0;
vi_ue_matrix = 0;
ve_ui_matrix = 0;
ve_ue_matrix = 0;
fe_face_values.reinit(cell, face_no);
fe_neighbor_face_values.reinit(neighbor, neighbor_face_no);
//double h = std::min(cell->diameter(),
// neighbor->diameter());
const double extent1 = cell-> measure() / cell-> face(face_no)-> measure();
const double extent2 = neighbor-> measure() / neighbor -> face(neighbor_face_no)-> measure();
const double penalty = compute_penalty(degree, extent1, extent2);
assemble_face_terms(fe_face_values,
fe_neighbor_face_values,
vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
typecase,
penalty);
//h);
neighbor->get_dof_indices(local_neighbor_dof_indices);
distribute_local_flux_to_global(vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
local_dof_indices,
local_neighbor_dof_indices);
}
}
}
}
constraints.distribute_local_to_global(local_matrix,
local_dof_indices,
system_matrix);
constraints.distribute_local_to_global(local_vector,
local_dof_indices,
system_rhs);
}
}
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
}
// i trial , j test
template<int dim>
void
FEMwDG<dim>::
assemble_cell_terms(
const FEValues<dim> &cell_fe,
FullMatrix<double> &cell_matrix,
Vector<double> &cell_vector,
const unsigned int &typecase,
const double &hk)
{
const unsigned int dofs_per_cell = cell_fe.dofs_per_cell;
const unsigned int n_q_points = cell_fe.n_quadrature_points;
// Information g and f
std::vector<double> g(n_q_points);
//std::vector<double> f(n_q_points);
boundary_function.value_list (cell_fe.get_quadrature_points(), g);
//rhs_function.value_list(cell_fe.get_quadrature_points(), f);
const FEValuesExtractors::Scalar verror(0);
const FEValuesExtractors::Scalar variable(1);
std::vector<double> rhs_values(n_q_points);
rhs_function.value_list(cell_fe.get_quadrature_points(),rhs_values);
for(unsigned int q=0; q<n_q_points; q++)
{
const Tensor<1,dim> beta_at_q_point = beta(cell_fe.quadrature_point(q));
const double kappa = coefficient<dim>(cell_fe.quadrature_point(q));
for(unsigned int i=0; i<dofs_per_cell; i++)
{
const Tensor <1, dim> grad_e_i = cell_fe[verror].gradient(i,q);
const Tensor <1, dim> grad_u_i = cell_fe[variable].gradient(i,q);
const double e_i = cell_fe[verror].value(i,q);
//const double phi_u_i = cell_fe[variable].value(i,q);
for(unsigned int j=0; j<dofs_per_cell; j++)
{
const Tensor <1, dim> grad_e_j = cell_fe[verror].gradient(j,q);
const Tensor <1, dim> grad_u_j = cell_fe[variable].gradient(j,q);
const double e_j = cell_fe[verror].value(j,q);
//const double phi_u_j = cell_fe[variable].value(j,q);
{
//---> contribution to B (SIP) - Domain
cell_matrix(i,j) += kappa * grad_u_j *
grad_e_i *
cell_fe.JxW(q);
//---> contribution to B (upw) - Domain
cell_matrix(i,j) += beta_at_q_point *
grad_u_j *
e_i *
cell_fe.JxW(q);
//---> contribution to BT (SIP) - Domain
cell_matrix(i,j) += kappa * grad_u_i *
grad_e_j *
cell_fe.JxW(q);
//---> contribution to BT (upw) - Domain
cell_matrix(i,j) += beta_at_q_point *
grad_u_i *
e_j *
cell_fe.JxW(q);
//---> contribution to G (upw) - Domain
// cell_matrix(i,j) += e_j *
// e_i *
// cell_fe.JxW(q); // <e,e1>
if (typecase == 1)
{
//<hk*(b*grad(e)),(b*grad(e1))> //---> ocultar cuando se quiera up -
cell_matrix(i,j)+= beta_at_q_point *
grad_e_j * hk *
grad_e_i *
beta_at_q_point *
cell_fe.JxW(q);
}
//---> contribution to G (SIP) - Domain
cell_matrix(i,j) += kappa * grad_e_j *
grad_e_i *
cell_fe.JxW(q);
}
}
cell_vector(i) += e_i *
rhs_values[q] *
cell_fe.JxW(q); // <f*e1>
}
}
}
template<int dim>
void
FEMwDG<dim>::
assemble_boundary_terms(
const FEFaceValues<dim> &face_fe,
FullMatrix<double> &local_matrix,
Vector<double> &local_vector,
const double &eta )
{
const unsigned int dofs_per_cell = face_fe.dofs_per_cell;
const unsigned int n_q_points = face_fe.n_quadrature_points;
const FEValuesExtractors::Scalar verror(0);
const FEValuesExtractors::Scalar variable(1);
// Information g and beta
std::vector<double> g(face_fe.n_quadrature_points);
boundary_function.value_list (face_fe.get_quadrature_points(), g);
const std::vector<Tensor<1,dim> > &normals = face_fe.get_normal_vectors ();
for(unsigned int q=0; q<n_q_points; q++)
{
const double kappa = coefficient<dim>(face_fe.quadrature_point(q));
const double beta_dot_n = beta(face_fe.quadrature_point(q)) * normals[q]; //
for(unsigned int i=0; i<dofs_per_cell; i++)
{
const Tensor <1, dim> grad_e_i = face_fe[verror].gradient(i,q);
const Tensor <1, dim> grad_u_i = face_fe[variable].gradient(i,q);
const double e_i = face_fe[verror].value(i,q);
const double u_i = face_fe[variable].value(i,q);
for(unsigned int j=0; j<dofs_per_cell; j++)
{
const Tensor <1, dim> grad_e_j = face_fe[verror].gradient(j,q);
const Tensor <1, dim> grad_u_j = face_fe[variable].gradient(j,q);
const double e_j = face_fe[verror].value(j,q);
const double u_j = face_fe[variable].value(j,q);
// //---> contribution to B (SIP) - Boundary
local_matrix(i,j) += ((kappa * eta) *
u_j *
e_i
- kappa *
grad_u_j * normals[q] *
e_i
- kappa *
u_j *
grad_e_i * normals[q]) *
face_fe.JxW(q);
//---> contribution to B (upw) - Boundary
local_matrix(i,j) += 0.5* (std::abs(beta_dot_n)-beta_dot_n)*
u_j *
e_i *
face_fe.JxW(q);
//---> contribution to BT (upw) - Boundary
local_matrix(i,j) += 0.5* (std::abs(beta_dot_n)-beta_dot_n)*
e_j *
u_i *
face_fe.JxW(q);
// //---> contribution to BT (SIP) - Boundary
local_matrix(i,j) += ((kappa * eta) *
e_j *
u_i
- kappa *
grad_e_j * normals[q] *
u_i
- kappa *
e_j *
grad_u_i * normals[q]) *
face_fe.JxW(q);
//---> contribution to G (upw) - Boundary
local_matrix(i,j) += 0.5*std::abs(beta_dot_n) *
e_j *
e_i *
face_fe.JxW(q);
//---> contribution to G (SIP) - Boundary
local_matrix(i,j) += (kappa * eta) *
e_j *
e_i *
face_fe.JxW(q);
}
//-----Vector L (SIP) (contribution - faces)--------
local_vector(i) += ((kappa * eta)* g[q] *
e_i
- g[q] *
kappa * grad_e_i * normals[q]) *
face_fe.JxW(q);
//------Vector L (upw) (contribution - faces)------
local_vector(i) += 0.5 *(std::abs(beta_dot_n)-beta_dot_n)*
g[q] *
e_i *
face_fe.JxW(q);
}
}
}
template<int dim>
void
FEMwDG<dim>::
assemble_face_terms(
const FEFaceValuesBase<dim> &fe_face_values,
const FEFaceValuesBase<dim> &fe_neighbor_face_values,
FullMatrix<double> &ui_vi_matrix,
FullMatrix<double> &ui_ve_matrix,
FullMatrix<double> &ue_vi_matrix,
FullMatrix<double> &ue_ve_matrix,
const unsigned int & typecase,
const double & eta)
{
const unsigned int n_face_points = fe_face_values.n_quadrature_points;
const unsigned int dofs_this_cell = fe_face_values.dofs_per_cell;
const unsigned int dofs_neighbor_cell = fe_neighbor_face_values.dofs_per_cell;
const std::vector<Tensor<1,dim> > &normals = fe_face_values.get_normal_vectors ();
//const std::vector<double> & JxW = fe_face_values.get_JxW_values();
const double penalty = (typecase==1) ? 1 : 0;
const FEValuesExtractors::Scalar verror(0);
const FEValuesExtractors::Scalar variable(1);
for(unsigned int q=0; q<n_face_points; q++)
{ //beta = beta dot n+
const double kappa = coefficient<dim>(fe_face_values.quadrature_point(q));
const double beta_dot_n = beta(fe_face_values.quadrature_point(q)) * normals[q]; //
for(unsigned int i=0; i<dofs_this_cell; i++)
{
const Tensor <1, dim> grad_ui_i = fe_face_values[variable].gradient(i,q);
const Tensor <1, dim> grad_ei_i = fe_face_values[verror].gradient(i,q);
const double ui_i = fe_face_values[variable].value(i,q);
const double ei_i = fe_face_values[verror].value(i,q);
for(unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor <1, dim> grad_ui_j = fe_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ei_j = fe_face_values[verror].gradient(j,q);
const double ui_j = fe_face_values[variable].value(j,q);
const double ei_j = fe_face_values[verror].value(j,q);
//---> contribution to B (SIP)
ui_vi_matrix(i,j) += (
(kappa* eta) *
ui_j *
ei_i
- (kappa*0.5)*
grad_ui_j * normals[q] *
ei_i
- (kappa*0.5)*
ui_j *
grad_ei_i * normals[q]
)*
fe_face_values.JxW(q);
// //---> contribution to B (upw)
ui_vi_matrix(i,j) += (-0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ui_j*
ei_i*
fe_face_values.JxW(q);
// //---> contribution to BT (upw)
ui_vi_matrix(i,j) += (-0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ei_j*
ui_i*
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ui_vi_matrix(i,j) += (
(kappa* eta) *
ei_j *
ui_i
- (kappa*0.5)*
grad_ei_j * normals[q] *
ui_i
- (kappa*0.5)*
ei_j *
grad_ui_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ui_vi_matrix(i,j) += penalty*0.5*std::abs(beta_dot_n)*
ei_j *
ei_i *
fe_face_values.JxW(q);
}
for(unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor <1, dim> grad_ue_j = fe_neighbor_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ee_j = fe_neighbor_face_values[verror].gradient(j,q);
const double ue_j = fe_neighbor_face_values[variable].value(j,q);
const double ee_j = fe_neighbor_face_values[verror].value(j,q);
//---> contribution to B (SIP)
ui_ve_matrix(i,j) += (-(kappa* eta)*
ue_j *
ei_i
- (kappa*0.5)*
grad_ue_j * normals[q] *
ei_i
+ (kappa*0.5)*
ue_j *
grad_ei_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to B (upw)
ui_ve_matrix(i,j) += (0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ue_j*
ei_i*
fe_face_values.JxW(q);
//---> contribution to BT (upw)
ui_ve_matrix(i,j) += (-0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ee_j*
ui_i*
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ui_ve_matrix(i,j) += (-(kappa* eta)*
ee_j *
ui_i
+ (kappa*0.5)*
grad_ee_j * normals[q] *
ui_i
- (kappa*0.5)*
ee_j *
grad_ui_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ui_ve_matrix(i,j) += -penalty*0.5*std::abs(beta_dot_n)*
ee_j *
ei_i*
fe_face_values.JxW(q);
//---> contribution to G (SIP)
ui_ve_matrix(i,j) += -(kappa * eta)*
ee_j *
ei_i *
fe_face_values.JxW(q);
}
}
for(unsigned int i=0; i<dofs_neighbor_cell; i++)
{
const Tensor <1, dim> grad_ue_i = fe_neighbor_face_values[variable].gradient(i,q);
const Tensor <1, dim> grad_ee_i = fe_neighbor_face_values[verror].gradient(i,q);
const double ue_i = fe_neighbor_face_values[variable].value(i,q);
const double ee_i = fe_neighbor_face_values[verror].value(i,q);
for(unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor <1, dim> grad_ui_j = fe_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ei_j = fe_face_values[verror].gradient(j,q);
const double ui_j = fe_face_values[variable].value(j,q);
const double ei_j = fe_face_values[verror].value(j,q);
ue_vi_matrix(i,j) += (- (kappa * eta)*
ui_j *
ee_i
+ (kappa*0.5)*
grad_ui_j * normals[q] *
ee_i
- (kappa*0.5)*
ui_j *
grad_ee_i * normals[q]
)*
fe_face_values.JxW(q);
//falta parte de BT sip
//---> contribution to B (upw)
ue_vi_matrix(i,j) += (-0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ui_j *
ee_i *
fe_face_values.JxW(q);
//---> contribution to BT (upw)
ue_vi_matrix(i,j) += (0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ei_j *
ue_i *
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ue_vi_matrix(i,j) += (- (kappa * eta)*
ei_j *
ue_i
- (kappa*0.5)*
grad_ei_j * normals[q] *
ue_i
+ (kappa*0.5)*
ei_j *
grad_ue_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ue_vi_matrix(i,j) += -penalty*0.5*std::abs(beta_dot_n)*
ei_j*
ee_i*
fe_face_values.JxW(q);
//---> contribution to G (SIP)
ue_vi_matrix(i,j) += -(kappa * eta)*
ei_j *
ee_i *
fe_face_values.JxW(q);
}
for(unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor <1, dim> grad_ue_j = fe_neighbor_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ee_j = fe_neighbor_face_values[verror].gradient(j,q);
const double ue_j = fe_neighbor_face_values[variable].value(j,q);
const double ee_j = fe_neighbor_face_values[verror].value(j,q);
// //---> contribution to B (SIP)
ue_ve_matrix(i,j) += ((kappa* eta)*
ue_j *
ee_i
+ (kappa*0.5)*
grad_ue_j * normals[q] *
ee_i
+ (kappa*0.5)*
ue_j *
grad_ee_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to B (upw)
ue_ve_matrix(i,j) += (0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ue_j *
ee_i*
fe_face_values.JxW(q);
//---> contribution to BT (upw)
ue_ve_matrix(i,j) += (0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ee_j *
ue_i*
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ue_ve_matrix(i,j) += ((kappa*eta)*
ee_j *
ue_i
+ (kappa*0.5)*
grad_ee_j * normals[q] *
ue_i
+ (kappa*0.5)*
ee_j *
grad_ue_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ue_ve_matrix(i,j) += penalty*0.5*std::abs(beta_dot_n)*
ee_j *
ee_i*
fe_face_values.JxW(q);
//---> contribution to G (SIP)
ue_ve_matrix(i,j) += (kappa* eta)*
ee_j *
ee_i *
fe_face_values.JxW(q);
}
}
}
}
template<int dim>
void
FEMwDG<dim>::
distribute_local_flux_to_global(
const FullMatrix<double> & vi_ui_matrix,
const FullMatrix<double> & vi_ue_matrix,
const FullMatrix<double> & ve_ui_matrix,
const FullMatrix<double> & ve_ue_matrix,
const std::vector<types::global_dof_index> & local_dof_indices,
const std::vector<types::global_dof_index> & local_neighbor_dof_indices)
{
constraints.distribute_local_to_global(vi_ui_matrix,
local_dof_indices,
system_matrix);
constraints.distribute_local_to_global(vi_ue_matrix,
local_dof_indices,
local_neighbor_dof_indices,
system_matrix);
constraints.distribute_local_to_global(ve_ui_matrix,
local_neighbor_dof_indices,
local_dof_indices,
system_matrix);
constraints.distribute_local_to_global(ve_ue_matrix,
local_neighbor_dof_indices,
system_matrix);
}
template <int dim>
void FEMwDG<dim>::
solve()
{
TimerOutput::Scope t(computing_timer, "solve");
LA::MPI::PreconditionJacobi prec_M;
{
LA::MPI::PreconditionJacobi::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
#endif
prec_M.initialize(system_matrix.block(0, 0), data);
}
TrilinosWrappers::PreconditionIdentity pre_aS;
{
TrilinosWrappers::PreconditionIdentity::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
#endif
pre_aS.initialize(system_matrix.block(0, 0), data);
}
auto &E = solution.block(0);
auto &U = solution.block(1);
const auto &L = system_rhs.block(0);
const auto M = linear_operator<LA::MPI::Vector>(system_matrix.block(0, 0));
const auto op_M = linear_operator(M, prec_M);
const auto op_B = linear_operator<LA::MPI::Vector>(system_matrix.block(0, 1));
const auto op_BT = linear_operator<LA::MPI::Vector>(system_matrix.block(1, 0));
const auto op_C = linear_operator<LA::MPI::Vector>(system_matrix.block(1, 1));
ReductionControl inner_solver_control2(5000,
1e-18 * system_rhs.l2_norm(),
1.e-2);
SolverCG<LA::MPI::Vector> cg(inner_solver_control2);
const auto op_M_inv = inverse_operator(M, cg, prec_M);
//const auto op_S = schur_complement(op_M_inv, op_B, op_BT, op_C);
const auto op_pre = linear_operator<LA::MPI::Vector>(prec_M);
const auto op_S = op_BT * op_M_inv * op_B;
const auto op_aS = op_BT * op_pre * op_B;
IterationNumberControl iteration_number_control_aS(30, 1.e-18);
SolverCG<LA::MPI::Vector> solver_aS(iteration_number_control_aS);
const auto preconditioner_S =
inverse_operator(op_S, solver_aS, pre_aS);
//inverse_operator(op_aS, solver_aS, pre_aS);
const auto schur_rhs = op_BT * op_M_inv * L ;
SolverControl solver_control_S(2000, 1.e-12);
SolverCG<LA::MPI::Vector> solver_S(solver_control_S);
const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );
// U = op_S_inv * schur_rhs;
// std::cout << solver_control_S.last_step()
// << " CG Schur complement iterations to obtain convergence."
// << std::endl;
// E = op_M_inv * (L - op_B * U);
}
template<int dim>
void
FEMwDG<dim>::
output_results(const unsigned int cycle) const
{
std::vector<std::string> solution_names;
solution_names = {"e","u"};
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(locally_relevant_solution,
solution_names);
Vector<float> subdomain(triangulation.n_active_cells());
for(unsigned int i=0; i<subdomain.size(); i++)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector(subdomain,"subdomain");
data_out.build_patches();
const std::string filename = ("solution." +
Utilities::int_to_string(
triangulation.locally_owned_subdomain(),4));
deallog << "Writing solution to <" << filename << ">" << std::endl;
std::ofstream output((filename + ".vtu").c_str());
data_out.write_vtu(output);
if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0 )
{
std::vector<std::string> filenames;
for(unsigned int i=0;
i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
i++)
{
filenames.push_back("solution." +
Utilities::int_to_string(i,4) +
".vtu");
}
std::ofstream master_output("solution.pvtu");
data_out.write_pvtu_record(master_output, filenames);
}
}
//-----------------------------------------EXACT SOLUTION--------------------------
template <int dim>
void FEMwDG<dim>::exact_plot(const unsigned int cycle)
{
VectorTools::project (dof_handler_u, constraints, QGauss<dim>(2),
ExactSolution<dim>(),
phi_0);
DataOut<dim> data_out;
const std::string filename2 = "exa-" + std::to_string(cycle)+ ".vtk";;
std::ofstream vtk_output (filename2.c_str());
data_out.attach_dof_handler(dof_handler_u);
data_out.add_data_vector(phi_0, "u_exact");
data_out.build_patches();
data_out.write_vtk(vtk_output);
}
//-----------------------------------------RUN CLASS--------------------------
template <int dim>
void FEMwDG<dim>::
run(const unsigned int refinement , const unsigned int typecase , const unsigned int n_cycles, const unsigned int refine_type)
{
Timer timer;
const std::string type_element1 =fe.get_name();
const std::string varDG_CG = "FESystem<" + std::to_string(dim) +">[FE_DGQ<" +std::to_string(dim)+">("+ std::to_string(degree)+
")-FE_Q<"+std::to_string(dim)+">("+ std::to_string(degree) +")]";
const std::string varDG_DG = "FESystem<" + std::to_string(dim) +">[FE_DGQ<" +std::to_string(dim)+">("+ std::to_string(degree)+
")-FE_DGQ<"+std::to_string(dim)+">("+ std::to_string(degree) +")]";
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
// cout << endl<< "Cycle number " <<cycle<<endl
// << "===========================================" << endl<< endl;
deallog << "Cycle: " << cycle << std::endl;
timer.start ();
make_grid(refinement,cycle, refine_type);
setup_system ();
// exact_plot(cycle);
Timer assemble_timer;
assemble_system (typecase);
std::cout << "Time of assemble_system: " << assemble_timer.cpu_time() << " seconds."
<< std::endl;
solve ();
// deallog << "CPUtime: " << timer.wall_time() << " seconds."<<std::endl;
// cout << "CPU time: " << timer.wall_time() << " seconds."<<std::endl;
// //estimate(typecase);
// output_results (cycle);
//compute_errors(cycle);
}
//tablas();
}
//-----------------------------------------MAIN--------------------------
int main(int argc, char *argv[])
{
try {
using namespace dealii;
const unsigned int dim = 2;
const unsigned int typecase = 1; //Up+=1 ; CentralFluxes=2;
const std::string filename = (typecase==1) ? "deallogUP" : "deallogCF";
std::ofstream logfile(filename);
deallog.attach(logfile);
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
numbers::invalid_unsigned_int);
const unsigned int refinement_initial = 6;
const unsigned int refine_type = 0; // ref. global = 0 ; ref. isotropico = 1 ;
const unsigned int Ptrial= 1; //Grado del polinomio (Ptrial = Ptest)
const unsigned int n_cycles=1;
FEMwDG<dim> Advection(Ptrial);
Advection.run(refinement_initial, typecase ,n_cycles, refine_type);
}
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;
}
return 0;
}
//ASFEM method by Juan Felipe Giraldo
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/generic_linear_algebra.h>
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
#include <deal.II/lac/block_linear_operator.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/linear_operator_tools.h>
#include <deal.II/lac/packaged_operation.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include<deal.II/lac/schur_complement.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.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_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <fstream>
#include <iostream>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.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/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/assembler.h>
#include <deal.II/meshworker/loop.h>
//#include "Functions2.cc"
#include <deal.II/base/function.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/lac/vector.h>
#include <cmath>
using namespace std;
using namespace dealii;
// ---------------------Beta ------------------------------
template <int dim>
Tensor<1,dim> beta (const Point<dim>)
{
Point<dim> field;
field(0) = 1;
field(1) = 0;
return field;
}
// ---------------------Kappa ------------------------------
template <int dim>
double coefficient(const Point<dim> &p)
{
return 1E-1;
}
template <int dim>
class RHS : public Function<dim>
{
public:
RHS() : Function<dim>(1)
{}
virtual double value(const Point<dim> &p,
const unsigned int component = 0 ) const;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues() : Function<dim>(1)
{}
virtual double value(const Point<dim> &p,
const unsigned int component = 0 ) const;
};
//----------------------------RHS -----------------------------------
template <int dim>
double
RHS<dim>::
value(const Point<dim> &p,
const unsigned int ) const
{
return 0.;
}
//-------- ----------Drichlet Boundary-------------------------------
template <int dim>
double BoundaryValues<dim>:: value(const Point<dim> &p,
const unsigned int ) const
{
const double x = p[0];
const double y = p[1];
//Eriksson
if (x < 1E-15)
{
const double kappa = coefficient<dim>(p);
double r1 = (1. + sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double r2 = (1. - sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
return (exp(r1*(x-1.)) - exp(r2*(x-1.)))/(exp(-r1)-exp(-r2)) * sin(M_PI*y);
}
else
{
return 0;
}
}
// -------------Exact Solution-------------------------------
template <int dim>
class ExactSolution : public Function<dim>
{
public:
ExactSolution() : Function<dim>() //first for the error, second for the scalar u value
{}
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override;
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
void ExactSolution<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
const double x = p[0];
const double y = p[1];
values(0) = 0; //error "exact solution"
const double kappa = coefficient<dim>(p);
double r1 = (1. + sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double r2 = (1. - sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double val = (exp(r1*(x-1.)) - exp(r2*(x-1.)))/(exp(-r1)-exp(-r2)) * sin(M_PI*y);
values(1) = val;
}
template <int dim>
double ExactSolution<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
const double x = p[0];
const double y = p[1];
const double kappa = coefficient<dim>(p);
double r1 = (1. + sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double r2 = (1. - sqrt(1.+4.*kappa*kappa*M_PI*M_PI))/(2.* kappa);
double val = (exp(r1*(x-1.)) - exp(r2*(x-1.)))/(exp(-r1)-exp(-r2)) * sin(M_PI*y);
return val;
}
//-----------------------------------CLASS TEMPLATE--FEMwDG-------------------------------------------
template <int dim>
class FEMwDG
{
public:
FEMwDG (const unsigned int degree);
void run (const unsigned int refinement_initial, const unsigned int typecase , const unsigned int n_cycles, const unsigned int refine_type);
~FEMwDG();
void run();
private:
using VectorType = TrilinosWrappers::MPI::Vector;
using MatrixType = TrilinosWrappers::SparseMatrix;
void make_grid(const unsigned int refinment_initial ,
const unsigned int cycle ,
const unsigned int refine_type);
void setup_system();
void exact_plot(const unsigned int cycle);
void assemble_system(const unsigned int typecase);
void assemble_cell_terms(const FEValues<dim> &cell_fe,
FullMatrix<double> &cell_matrix,
Vector<double> &cell_vector,
const unsigned int &typecase,
const double & h);
void assemble_boundary_terms(const FEFaceValues<dim> &face_fe,
FullMatrix<double> &local_matrix,
Vector<double> &local_vector,
const double &h);
//const unsigned int &typecase,
void assemble_face_terms(const FEFaceValuesBase<dim> &fe_face_values,
const FEFaceValuesBase<dim> &fe_neighbor_face_values,
FullMatrix<double> &vi_ui_matrix,
FullMatrix<double> &vi_ue_matrix,
FullMatrix<double> &ve_ui_matrix,
FullMatrix<double> &ve_ue_matrix,
const unsigned int &typecase,
const double & h);
void estimate(const unsigned int typecase);
void distribute_local_flux_to_global(
const FullMatrix<double> & vi_ui_matrix,
const FullMatrix<double> & vi_ue_matrix,
const FullMatrix<double> & ve_ui_matrix,
const FullMatrix<double> & ve_ue_matrix,
const std::vector<types::global_dof_index> & local_dof_indices,
const std::vector<types::global_dof_index> & local_neighbor_dof_indices);
void solve();
//void estimate(const unsigned int typecase);
void output_results(const unsigned int cycle) const;
const unsigned int degree;
//const unsignet int typecase;
parallel::distributed::Triangulation<dim> triangulation;
FESystem<dim> fe;
FE_DGQ<dim> fe_e;
FE_Q<dim> fe_u;
DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_e;
DoFHandler<dim> dof_handler_u;
LA::MPI::BlockSparseMatrix system_matrix;
LA::MPI::BlockSparseMatrix preconditioner_matrix;
LA::MPI::BlockVector locally_relevant_solution;
LA::MPI::BlockVector system_rhs;
LA::MPI::BlockVector solution;
std::vector<IndexSet> owned_partitioning;
std::vector<IndexSet> relevant_partitioning;
AffineConstraints<double> constraints;
AffineConstraints<double> constraints_u;
Vector<double> estimated_error_per_cell;
//BlockVector<double> estimates;
Vector<double> phi_0;
ConditionalOStream pcout;
TimerOutput computing_timer;
SolverControl solver_control;
TrilinosWrappers::SolverDirect solver;
const RHS<dim> rhs_function;
const BoundaryValues<dim> boundary_function;
const ExactSolution<dim> exact_solution;
};
template <int dim>
FEMwDG<dim>::
FEMwDG(const unsigned int degree)
:
degree(degree),
//n_refine(n_refine),
triangulation(MPI_COMM_WORLD,
typename Triangulation<dim>::MeshSmoothing
(Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening)),
fe(FE_DGQ<dim>(degree), 1,FE_Q<dim>(degree), 1), // test(e)->(D.G)- trial(u)->(C.G). Second term refers to a number of copies of fe, copy for direction (in this case is 1).
//fe(FE_DGQ<dim>(degree), 1,FE_DGQ<dim>(degree), 1), // --> ONLY DG
fe_e(degree),
fe_u(degree),
//fe(FESystem<dim>(FE_DGQ<dim>(degree), dim) , 1,FE_DGQ<dim>(degree),1),
dof_handler(triangulation),
dof_handler_e (triangulation),
dof_handler_u (triangulation),
pcout(std::cout,
Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
computing_timer(MPI_COMM_WORLD,
pcout,
TimerOutput::summary,
TimerOutput::wall_times),
solver_control(1),
solver(solver_control),
rhs_function(),
boundary_function()
//estimates(1)
{}
template <int dim>
FEMwDG<dim>::
~FEMwDG()
{
dof_handler.clear();
}
//-----------------------------------------MESH CREATION --------------------------
template <int dim>
void FEMwDG<dim>::make_grid(const unsigned int refinement_initial , const unsigned int cycle, const unsigned int refine_type)
{
if (cycle==0)
{
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(refinement_initial);
}
else
{
if (refine_type==0) // Global refinement
triangulation.refine_global(1);
if (refine_type==1) // Adaptive refinement (Dorfler marking)
{triangulation.refine_global(1);
// parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
//triangulation, estimated_error_per_cell, 0.3, 0.1);
//GridRefinement::refine_and_coarsen_fixed_fraction(
//triangulation, estimates.block(0), 0.5, 0.3);
//triangulation.execute_coarsening_and_refinement();
}
}
}
//-----------------------------------------RUN CLASS--------------------------
template <int dim>
void
FEMwDG<dim>::
setup_system()
{
TimerOutput::Scope t(computing_timer, "setup");
dof_handler.distribute_dofs(fe);
DoFRenumbering::component_wise(dof_handler);
std::vector<types::global_dof_index> dofs_per_block(2); //Grados de libertad por componentes. (escalar - escalar)
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block);//, block_component);
const unsigned int n_e = dofs_per_block[0],
n_u = dofs_per_block[1];
//std::cout << "Number of active cells: " << triangulation.n_active_cells()
// << std::endl
std::cout << "Total number of cells: " << triangulation.n_active_cells()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << n_e << '+' << n_u << ')' << std::endl;
// set the configuration for the projection - exact solution
dof_handler_u.distribute_dofs(fe_u);
DynamicSparsityPattern dsp_u (n_u, n_u);
phi_0.reinit (n_u);
constraints_u.clear ();
DoFTools::make_flux_sparsity_pattern (dof_handler_u, dsp_u, constraints_u, false);
DoFRenumbering::component_wise (dof_handler_u);
owned_partitioning.resize(2);
owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, n_e);
owned_partitioning[1] = dof_handler.locally_owned_dofs().get_view(n_e, n_e + n_u);
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
relevant_partitioning.resize(2);
relevant_partitioning[0] = locally_relevant_dofs.get_view(0, n_e);
relevant_partitioning[1] = locally_relevant_dofs.get_view(n_e, n_e + n_u);
constraints.clear();
constraints.close();
system_matrix.clear();
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
DoFTools::make_flux_sparsity_pattern(dof_handler,
dsp);
SparsityTools::distribute_sparsity_pattern(
dsp,
dof_handler.locally_owned_dofs(),
MPI_COMM_WORLD,
locally_relevant_dofs);
system_matrix.reinit(owned_partitioning,
dsp,
MPI_COMM_WORLD);
preconditioner_matrix.clear();
DoFTools::make_flux_sparsity_pattern(dof_handler,
dsp);
SparsityTools::distribute_sparsity_pattern(
dsp,
Utilities::MPI::all_gather(MPI_COMM_WORLD,dof_handler.locally_owned_dofs()),
MPI_COMM_WORLD,
locally_relevant_dofs);
preconditioner_matrix.reinit(owned_partitioning,
dsp,
MPI_COMM_WORLD);
locally_relevant_solution.reinit(owned_partitioning,
relevant_partitioning,
MPI_COMM_WORLD);
system_rhs.reinit(owned_partitioning, MPI_COMM_WORLD);
solution.reinit(owned_partitioning, MPI_COMM_WORLD);
}
double compute_penalty(const unsigned int fe_degree,
const double cell_extend_left,
const double cell_extend_right)
{
const double degree = std::max(1., static_cast<double>(fe_degree));
return degree * (degree + 1.) * 0.5 *
(1. / cell_extend_left + 1. / cell_extend_right);
}
template <int dim>
void
FEMwDG<dim>::
assemble_system(const unsigned int typecase)
{
TimerOutput::Scope t(computing_timer, "assembly");
QGauss<dim> quadrature_formula(fe.degree+1);
QGauss<dim-1> face_quadrature_formula(fe.degree+1);
const UpdateFlags update_flags = update_values
| update_gradients
| update_quadrature_points
| update_JxW_values;
const UpdateFlags face_update_flags = update_values
| update_gradients
| update_normal_vectors
| update_quadrature_points
| update_JxW_values;
const unsigned int dofs_per_cell = fe.dofs_per_cell;
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<types::global_dof_index> local_neighbor_dof_indices(dofs_per_cell);
FEValues<dim> fe_values(fe, quadrature_formula, update_flags);
FEFaceValues<dim> fe_face_values(fe,face_quadrature_formula,
face_update_flags);
FEFaceValues<dim> fe_neighbor_face_values(fe, face_quadrature_formula,
face_update_flags);
FESubfaceValues<dim> fe_subface_values(fe, face_quadrature_formula,
face_update_flags);
FullMatrix<double> local_matrix(dofs_per_cell,dofs_per_cell);
Vector<double> local_vector(dofs_per_cell);
FullMatrix<double> vi_ui_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> vi_ue_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> ve_ui_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> ve_ue_matrix(dofs_per_cell, dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for(; cell!=endc; cell++)
{
if(cell->is_locally_owned())
{
local_matrix = 0;
local_vector = 0;
fe_values.reinit(cell);
double h = cell->diameter();
//------------DOMAIN INTEGRATOR----------------------------------------------
assemble_cell_terms(fe_values,
local_matrix,
local_vector,
typecase,
h);
cell->get_dof_indices(local_dof_indices);
for(unsigned int face_no=0;
face_no< GeometryInfo<dim>::faces_per_cell;
face_no++)
{
//-------------EXTERNAL BOUNDARY INTEGRATOR-------------------------------------
typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
if(face->at_boundary() )
{
const double extent1 = cell->measure() / cell->face(face_no)->measure();
const double penalty = compute_penalty(degree, extent1, extent1);
fe_face_values.reinit(cell, face_no);
assemble_boundary_terms(fe_face_values,
local_matrix,
local_vector,
penalty);
//typecase,h);
}
//-------------SKELETON INTEGRATOR ----------------------------------------------
else
{
Assert(cell->neighbor(face_no).state() == IteratorState::valid,
ExcInternalError());
typename DoFHandler<dim>::cell_iterator neighbor =
cell->neighbor(face_no);
if(face->has_children())
{
const unsigned int neighbor_face_no =
cell->neighbor_of_neighbor(face_no);
for(unsigned int subface_no=0;
subface_no < face->number_of_children();
++subface_no)
{
typename DoFHandler<dim>::cell_iterator neighbor_child =
cell->neighbor_child_on_subface(face_no, subface_no);
Assert(!neighbor_child->has_children(), ExcInternalError());
vi_ui_matrix = 0;
vi_ue_matrix = 0;
ve_ui_matrix = 0;
ve_ue_matrix = 0;
fe_subface_values.reinit(cell, face_no, subface_no);
fe_neighbor_face_values.reinit(neighbor_child,
neighbor_face_no);
const double extent1 = cell-> measure() / cell-> face(face_no)-> measure();
const double extent2 = neighbor_child -> measure() / neighbor_child -> face(subface_no)-> measure();
const double penalty = compute_penalty(degree, extent1, extent2);
assemble_face_terms(fe_subface_values,
fe_neighbor_face_values,
vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
typecase,
penalty);
//h);
neighbor_child->get_dof_indices(local_neighbor_dof_indices);
distribute_local_flux_to_global(vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
local_dof_indices,
local_neighbor_dof_indices);
}
}
else
{
//---------------ENSAMBLANDO CONTRIBUCIÓN LOCAL EN LOS BORDES AL SISTEMA GLOBAL CUANDO EXISTE VECINO-------------------------------------
if(neighbor->level() == cell->level() &&
neighbor->index() > cell->index() )
{
const unsigned int neighbor_face_no =
cell->neighbor_of_neighbor(face_no);
vi_ui_matrix = 0;
vi_ue_matrix = 0;
ve_ui_matrix = 0;
ve_ue_matrix = 0;
fe_face_values.reinit(cell, face_no);
fe_neighbor_face_values.reinit(neighbor, neighbor_face_no);
//double h = std::min(cell->diameter(),
// neighbor->diameter());
const double extent1 = cell-> measure() / cell-> face(face_no)-> measure();
const double extent2 = neighbor-> measure() / neighbor -> face(neighbor_face_no)-> measure();
const double penalty = compute_penalty(degree, extent1, extent2);
assemble_face_terms(fe_face_values,
fe_neighbor_face_values,
vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
typecase,
penalty);
//h);
neighbor->get_dof_indices(local_neighbor_dof_indices);
distribute_local_flux_to_global(vi_ui_matrix,
vi_ue_matrix,
ve_ui_matrix,
ve_ue_matrix,
local_dof_indices,
local_neighbor_dof_indices);
}
}
}
}
constraints.distribute_local_to_global(local_matrix,
local_dof_indices,
system_matrix);
constraints.distribute_local_to_global(local_vector,
local_dof_indices,
system_rhs);
}
}
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
}
// i trial , j test
template<int dim>
void
FEMwDG<dim>::
assemble_cell_terms(
const FEValues<dim> &cell_fe,
FullMatrix<double> &cell_matrix,
Vector<double> &cell_vector,
const unsigned int &typecase,
const double &hk)
{
const unsigned int dofs_per_cell = cell_fe.dofs_per_cell;
const unsigned int n_q_points = cell_fe.n_quadrature_points;
// Information g and f
std::vector<double> g(n_q_points);
//std::vector<double> f(n_q_points);
boundary_function.value_list (cell_fe.get_quadrature_points(), g);
//rhs_function.value_list(cell_fe.get_quadrature_points(), f);
const FEValuesExtractors::Scalar verror(0);
const FEValuesExtractors::Scalar variable(1);
std::vector<double> rhs_values(n_q_points);
rhs_function.value_list(cell_fe.get_quadrature_points(),rhs_values);
for(unsigned int q=0; q<n_q_points; q++)
{
const Tensor<1,dim> beta_at_q_point = beta(cell_fe.quadrature_point(q));
const double kappa = coefficient<dim>(cell_fe.quadrature_point(q));
for(unsigned int i=0; i<dofs_per_cell; i++)
{
const Tensor <1, dim> grad_e_i = cell_fe[verror].gradient(i,q);
const Tensor <1, dim> grad_u_i = cell_fe[variable].gradient(i,q);
const double e_i = cell_fe[verror].value(i,q);
//const double phi_u_i = cell_fe[variable].value(i,q);
for(unsigned int j=0; j<dofs_per_cell; j++)
{
const Tensor <1, dim> grad_e_j = cell_fe[verror].gradient(j,q);
const Tensor <1, dim> grad_u_j = cell_fe[variable].gradient(j,q);
const double e_j = cell_fe[verror].value(j,q);
//const double phi_u_j = cell_fe[variable].value(j,q);
{
//---> contribution to B (SIP) - Domain
cell_matrix(i,j) += kappa * grad_u_j *
grad_e_i *
cell_fe.JxW(q);
//---> contribution to B (upw) - Domain
cell_matrix(i,j) += beta_at_q_point *
grad_u_j *
e_i *
cell_fe.JxW(q);
//---> contribution to BT (SIP) - Domain
cell_matrix(i,j) += kappa * grad_u_i *
grad_e_j *
cell_fe.JxW(q);
//---> contribution to BT (upw) - Domain
cell_matrix(i,j) += beta_at_q_point *
grad_u_i *
e_j *
cell_fe.JxW(q);
//---> contribution to G (upw) - Domain
// cell_matrix(i,j) += e_j *
// e_i *
// cell_fe.JxW(q); // <e,e1>
if (typecase == 1)
{
//<hk*(b*grad(e)),(b*grad(e1))> //---> ocultar cuando se quiera up -
cell_matrix(i,j)+= beta_at_q_point *
grad_e_j * hk *
grad_e_i *
beta_at_q_point *
cell_fe.JxW(q);
}
//---> contribution to G (SIP) - Domain
cell_matrix(i,j) += kappa * grad_e_j *
grad_e_i *
cell_fe.JxW(q);
}
}
cell_vector(i) += e_i *
rhs_values[q] *
cell_fe.JxW(q); // <f*e1>
}
}
}
template<int dim>
void
FEMwDG<dim>::
assemble_boundary_terms(
const FEFaceValues<dim> &face_fe,
FullMatrix<double> &local_matrix,
Vector<double> &local_vector,
const double &eta )
{
const unsigned int dofs_per_cell = face_fe.dofs_per_cell;
const unsigned int n_q_points = face_fe.n_quadrature_points;
const FEValuesExtractors::Scalar verror(0);
const FEValuesExtractors::Scalar variable(1);
// Information g and beta
std::vector<double> g(face_fe.n_quadrature_points);
boundary_function.value_list (face_fe.get_quadrature_points(), g);
const std::vector<Tensor<1,dim> > &normals = face_fe.get_normal_vectors ();
for(unsigned int q=0; q<n_q_points; q++)
{
const double kappa = coefficient<dim>(face_fe.quadrature_point(q));
const double beta_dot_n = beta(face_fe.quadrature_point(q)) * normals[q]; //
for(unsigned int i=0; i<dofs_per_cell; i++)
{
const Tensor <1, dim> grad_e_i = face_fe[verror].gradient(i,q);
const Tensor <1, dim> grad_u_i = face_fe[variable].gradient(i,q);
const double e_i = face_fe[verror].value(i,q);
const double u_i = face_fe[variable].value(i,q);
for(unsigned int j=0; j<dofs_per_cell; j++)
{
const Tensor <1, dim> grad_e_j = face_fe[verror].gradient(j,q);
const Tensor <1, dim> grad_u_j = face_fe[variable].gradient(j,q);
const double e_j = face_fe[verror].value(j,q);
const double u_j = face_fe[variable].value(j,q);
// //---> contribution to B (SIP) - Boundary
local_matrix(i,j) += ((kappa * eta) *
u_j *
e_i
- kappa *
grad_u_j * normals[q] *
e_i
- kappa *
u_j *
grad_e_i * normals[q]) *
face_fe.JxW(q);
//---> contribution to B (upw) - Boundary
local_matrix(i,j) += 0.5* (std::abs(beta_dot_n)-beta_dot_n)*
u_j *
e_i *
face_fe.JxW(q);
//---> contribution to BT (upw) - Boundary
local_matrix(i,j) += 0.5* (std::abs(beta_dot_n)-beta_dot_n)*
e_j *
u_i *
face_fe.JxW(q);
// //---> contribution to BT (SIP) - Boundary
local_matrix(i,j) += ((kappa * eta) *
e_j *
u_i
- kappa *
grad_e_j * normals[q] *
u_i
- kappa *
e_j *
grad_u_i * normals[q]) *
face_fe.JxW(q);
//---> contribution to G (upw) - Boundary
local_matrix(i,j) += 0.5*std::abs(beta_dot_n) *
e_j *
e_i *
face_fe.JxW(q);
//---> contribution to G (SIP) - Boundary
local_matrix(i,j) += (kappa * eta) *
e_j *
e_i *
face_fe.JxW(q);
}
//-----Vector L (SIP) (contribution - faces)--------
local_vector(i) += ((kappa * eta)* g[q] *
e_i
- g[q] *
kappa * grad_e_i * normals[q]) *
face_fe.JxW(q);
//------Vector L (upw) (contribution - faces)------
local_vector(i) += 0.5 *(std::abs(beta_dot_n)-beta_dot_n)*
g[q] *
e_i *
face_fe.JxW(q);
}
}
}
template<int dim>
void
FEMwDG<dim>::
assemble_face_terms(
const FEFaceValuesBase<dim> &fe_face_values,
const FEFaceValuesBase<dim> &fe_neighbor_face_values,
FullMatrix<double> &ui_vi_matrix,
FullMatrix<double> &ui_ve_matrix,
FullMatrix<double> &ue_vi_matrix,
FullMatrix<double> &ue_ve_matrix,
const unsigned int & typecase,
const double & eta)
{
const unsigned int n_face_points = fe_face_values.n_quadrature_points;
const unsigned int dofs_this_cell = fe_face_values.dofs_per_cell;
const unsigned int dofs_neighbor_cell = fe_neighbor_face_values.dofs_per_cell;
const std::vector<Tensor<1,dim> > &normals = fe_face_values.get_normal_vectors ();
//const std::vector<double> & JxW = fe_face_values.get_JxW_values();
const double penalty = (typecase==1) ? 1 : 0;
const FEValuesExtractors::Scalar verror(0);
const FEValuesExtractors::Scalar variable(1);
for(unsigned int q=0; q<n_face_points; q++)
{ //beta = beta dot n+
const double kappa = coefficient<dim>(fe_face_values.quadrature_point(q));
const double beta_dot_n = beta(fe_face_values.quadrature_point(q)) * normals[q]; //
for(unsigned int i=0; i<dofs_this_cell; i++)
{
const Tensor <1, dim> grad_ui_i = fe_face_values[variable].gradient(i,q);
const Tensor <1, dim> grad_ei_i = fe_face_values[verror].gradient(i,q);
const double ui_i = fe_face_values[variable].value(i,q);
const double ei_i = fe_face_values[verror].value(i,q);
for(unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor <1, dim> grad_ui_j = fe_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ei_j = fe_face_values[verror].gradient(j,q);
const double ui_j = fe_face_values[variable].value(j,q);
const double ei_j = fe_face_values[verror].value(j,q);
//---> contribution to B (SIP)
ui_vi_matrix(i,j) += (
(kappa* eta) *
ui_j *
ei_i
- (kappa*0.5)*
grad_ui_j * normals[q] *
ei_i
- (kappa*0.5)*
ui_j *
grad_ei_i * normals[q]
)*
fe_face_values.JxW(q);
// //---> contribution to B (upw)
ui_vi_matrix(i,j) += (-0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ui_j*
ei_i*
fe_face_values.JxW(q);
// //---> contribution to BT (upw)
ui_vi_matrix(i,j) += (-0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ei_j*
ui_i*
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ui_vi_matrix(i,j) += (
(kappa* eta) *
ei_j *
ui_i
- (kappa*0.5)*
grad_ei_j * normals[q] *
ui_i
- (kappa*0.5)*
ei_j *
grad_ui_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ui_vi_matrix(i,j) += penalty*0.5*std::abs(beta_dot_n)*
ei_j *
ei_i *
fe_face_values.JxW(q);
}
for(unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor <1, dim> grad_ue_j = fe_neighbor_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ee_j = fe_neighbor_face_values[verror].gradient(j,q);
const double ue_j = fe_neighbor_face_values[variable].value(j,q);
const double ee_j = fe_neighbor_face_values[verror].value(j,q);
//---> contribution to B (SIP)
ui_ve_matrix(i,j) += (-(kappa* eta)*
ue_j *
ei_i
- (kappa*0.5)*
grad_ue_j * normals[q] *
ei_i
+ (kappa*0.5)*
ue_j *
grad_ei_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to B (upw)
ui_ve_matrix(i,j) += (0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ue_j*
ei_i*
fe_face_values.JxW(q);
//---> contribution to BT (upw)
ui_ve_matrix(i,j) += (-0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ee_j*
ui_i*
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ui_ve_matrix(i,j) += (-(kappa* eta)*
ee_j *
ui_i
+ (kappa*0.5)*
grad_ee_j * normals[q] *
ui_i
- (kappa*0.5)*
ee_j *
grad_ui_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ui_ve_matrix(i,j) += -penalty*0.5*std::abs(beta_dot_n)*
ee_j *
ei_i*
fe_face_values.JxW(q);
//---> contribution to G (SIP)
ui_ve_matrix(i,j) += -(kappa * eta)*
ee_j *
ei_i *
fe_face_values.JxW(q);
}
}
for(unsigned int i=0; i<dofs_neighbor_cell; i++)
{
const Tensor <1, dim> grad_ue_i = fe_neighbor_face_values[variable].gradient(i,q);
const Tensor <1, dim> grad_ee_i = fe_neighbor_face_values[verror].gradient(i,q);
const double ue_i = fe_neighbor_face_values[variable].value(i,q);
const double ee_i = fe_neighbor_face_values[verror].value(i,q);
for(unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor <1, dim> grad_ui_j = fe_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ei_j = fe_face_values[verror].gradient(j,q);
const double ui_j = fe_face_values[variable].value(j,q);
const double ei_j = fe_face_values[verror].value(j,q);
ue_vi_matrix(i,j) += (- (kappa * eta)*
ui_j *
ee_i
+ (kappa*0.5)*
grad_ui_j * normals[q] *
ee_i
- (kappa*0.5)*
ui_j *
grad_ee_i * normals[q]
)*
fe_face_values.JxW(q);
//falta parte de BT sip
//---> contribution to B (upw)
ue_vi_matrix(i,j) += (-0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ui_j *
ee_i *
fe_face_values.JxW(q);
//---> contribution to BT (upw)
ue_vi_matrix(i,j) += (0.5*beta_dot_n - 0.5*penalty*std::abs(beta_dot_n))*
ei_j *
ue_i *
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ue_vi_matrix(i,j) += (- (kappa * eta)*
ei_j *
ue_i
- (kappa*0.5)*
grad_ei_j * normals[q] *
ue_i
+ (kappa*0.5)*
ei_j *
grad_ue_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ue_vi_matrix(i,j) += -penalty*0.5*std::abs(beta_dot_n)*
ei_j*
ee_i*
fe_face_values.JxW(q);
//---> contribution to G (SIP)
ue_vi_matrix(i,j) += -(kappa * eta)*
ei_j *
ee_i *
fe_face_values.JxW(q);
}
for(unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor <1, dim> grad_ue_j = fe_neighbor_face_values[variable].gradient(j,q);
const Tensor <1, dim> grad_ee_j = fe_neighbor_face_values[verror].gradient(j,q);
const double ue_j = fe_neighbor_face_values[variable].value(j,q);
const double ee_j = fe_neighbor_face_values[verror].value(j,q);
// //---> contribution to B (SIP)
ue_ve_matrix(i,j) += ((kappa* eta)*
ue_j *
ee_i
+ (kappa*0.5)*
grad_ue_j * normals[q] *
ee_i
+ (kappa*0.5)*
ue_j *
grad_ee_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to B (upw)
ue_ve_matrix(i,j) += (0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ue_j *
ee_i*
fe_face_values.JxW(q);
//---> contribution to BT (upw)
ue_ve_matrix(i,j) += (0.5*beta_dot_n + 0.5*penalty*std::abs(beta_dot_n))*
ee_j *
ue_i*
fe_face_values.JxW(q);
//---> contribution to BT (SIP)
ue_ve_matrix(i,j) += ((kappa*eta)*
ee_j *
ue_i
+ (kappa*0.5)*
grad_ee_j * normals[q] *
ue_i
+ (kappa*0.5)*
ee_j *
grad_ue_i * normals[q]
)*
fe_face_values.JxW(q);
//---> contribution to G (upw)
ue_ve_matrix(i,j) += penalty*0.5*std::abs(beta_dot_n)*
ee_j *
ee_i*
fe_face_values.JxW(q);
//---> contribution to G (SIP)
ue_ve_matrix(i,j) += (kappa* eta)*
ee_j *
ee_i *
fe_face_values.JxW(q);
}
}
}
}
template<int dim>
void
FEMwDG<dim>::
distribute_local_flux_to_global(
const FullMatrix<double> & vi_ui_matrix,
const FullMatrix<double> & vi_ue_matrix,
const FullMatrix<double> & ve_ui_matrix,
const FullMatrix<double> & ve_ue_matrix,
const std::vector<types::global_dof_index> & local_dof_indices,
const std::vector<types::global_dof_index> & local_neighbor_dof_indices)
{
constraints.distribute_local_to_global(vi_ui_matrix,
local_dof_indices,
system_matrix);
constraints.distribute_local_to_global(vi_ue_matrix,
local_dof_indices,
local_neighbor_dof_indices,
system_matrix);
constraints.distribute_local_to_global(ve_ui_matrix,
local_neighbor_dof_indices,
local_dof_indices,
system_matrix);
constraints.distribute_local_to_global(ve_ue_matrix,
local_neighbor_dof_indices,
system_matrix);
}
template <int dim>
void FEMwDG<dim>::
solve()
{
TimerOutput::Scope t(computing_timer, "solve");
LA::MPI::PreconditionJacobi prec_M;
{
LA::MPI::PreconditionJacobi::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
#endif
prec_M.initialize(system_matrix.block(0, 0), data);
}
TrilinosWrappers::PreconditionIdentity pre_aS;
{
TrilinosWrappers::PreconditionIdentity::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
#endif
pre_aS.initialize(system_matrix.block(0, 0), data);
}
auto &E = solution.block(0);
auto &U = solution.block(1);
const auto &L = system_rhs.block(0);
const auto M = linear_operator<LA::MPI::Vector>(system_matrix.block(0, 0));
const auto op_M = linear_operator(M, prec_M);
const auto op_B = linear_operator<LA::MPI::Vector>(system_matrix.block(0, 1));
const auto op_BT = linear_operator<LA::MPI::Vector>(system_matrix.block(1, 0));
const auto op_C = linear_operator<LA::MPI::Vector>(system_matrix.block(1, 1));
ReductionControl inner_solver_control2(5000,
1e-18 * system_rhs.l2_norm(),
1.e-2);
SolverCG<LA::MPI::Vector> cg(inner_solver_control2);
const auto op_M_inv = inverse_operator(M, cg, prec_M);
//const auto op_S = schur_complement(op_M_inv, op_B, op_BT, op_C);
const auto op_pre = linear_operator<LA::MPI::Vector>(prec_M);
const auto op_S = op_BT * op_M_inv * op_B;
const auto op_aS = op_BT * op_pre * op_B;
IterationNumberControl iteration_number_control_aS(30, 1.e-18);
SolverCG<LA::MPI::Vector> solver_aS(iteration_number_control_aS);
const auto preconditioner_S =
//inverse_operator(op_S, solver_aS, pre_aS);
inverse_operator(op_aS, solver_aS, pre_aS);
const auto schur_rhs = op_BT * op_M_inv * L ;
SolverControl solver_control_S(2000, 1.e-12);
SolverCG<LA::MPI::Vector> solver_S(solver_control_S);
const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );
// U = op_S_inv * schur_rhs;
// std::cout << solver_control_S.last_step()
// << " CG Schur complement iterations to obtain convergence."
// << std::endl;
// E = op_M_inv * (L - op_B * U);
}
template<int dim>
void
FEMwDG<dim>::
output_results(const unsigned int cycle) const
{
std::vector<std::string> solution_names;
solution_names = {"e","u"};
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(locally_relevant_solution,
solution_names);
Vector<float> subdomain(triangulation.n_active_cells());
for(unsigned int i=0; i<subdomain.size(); i++)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector(subdomain,"subdomain");
data_out.build_patches();
const std::string filename = ("solution." +
Utilities::int_to_string(
triangulation.locally_owned_subdomain(),4));
deallog << "Writing solution to <" << filename << ">" << std::endl;
std::ofstream output((filename + ".vtu").c_str());
data_out.write_vtu(output);
if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0 )
{
std::vector<std::string> filenames;
for(unsigned int i=0;
i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
i++)
{
filenames.push_back("solution." +
Utilities::int_to_string(i,4) +
".vtu");
}
std::ofstream master_output("solution.pvtu");
data_out.write_pvtu_record(master_output, filenames);
}
}
//-----------------------------------------EXACT SOLUTION--------------------------
template <int dim>
void FEMwDG<dim>::exact_plot(const unsigned int cycle)
{
VectorTools::project (dof_handler_u, constraints, QGauss<dim>(2),
ExactSolution<dim>(),
phi_0);
DataOut<dim> data_out;
const std::string filename2 = "exa-" + std::to_string(cycle)+ ".vtk";;
std::ofstream vtk_output (filename2.c_str());
data_out.attach_dof_handler(dof_handler_u);
data_out.add_data_vector(phi_0, "u_exact");
data_out.build_patches();
data_out.write_vtk(vtk_output);
}
//-----------------------------------------RUN CLASS--------------------------
template <int dim>
void FEMwDG<dim>::
run(const unsigned int refinement , const unsigned int typecase , const unsigned int n_cycles, const unsigned int refine_type)
{
Timer timer;
const std::string type_element1 =fe.get_name();
const std::string varDG_CG = "FESystem<" + std::to_string(dim) +">[FE_DGQ<" +std::to_string(dim)+">("+ std::to_string(degree)+
")-FE_Q<"+std::to_string(dim)+">("+ std::to_string(degree) +")]";
const std::string varDG_DG = "FESystem<" + std::to_string(dim) +">[FE_DGQ<" +std::to_string(dim)+">("+ std::to_string(degree)+
")-FE_DGQ<"+std::to_string(dim)+">("+ std::to_string(degree) +")]";
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
// cout << endl<< "Cycle number " <<cycle<<endl
// << "===========================================" << endl<< endl;
deallog << "Cycle: " << cycle << std::endl;
timer.start ();
make_grid(refinement,cycle, refine_type);
setup_system ();
// exact_plot(cycle);
Timer assemble_timer;
assemble_system (typecase);
std::cout << "Time of assemble_system: " << assemble_timer.cpu_time() << " seconds."
<< std::endl;
solve ();
// deallog << "CPUtime: " << timer.wall_time() << " seconds."<<std::endl;
// cout << "CPU time: " << timer.wall_time() << " seconds."<<std::endl;
// //estimate(typecase);
// output_results (cycle);
//compute_errors(cycle);
}
//tablas();
}
//-----------------------------------------MAIN--------------------------
int main(int argc, char *argv[])
{
try {
using namespace dealii;
const unsigned int dim = 2;
const unsigned int typecase = 1; //Up+=1 ; CentralFluxes=2;
const std::string filename = (typecase==1) ? "deallogUP" : "deallogCF";
std::ofstream logfile(filename);
deallog.attach(logfile);
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
numbers::invalid_unsigned_int);
const unsigned int refinement_initial = 6;
const unsigned int refine_type = 0; // ref. global = 0 ; ref. isotropico = 1 ;
const unsigned int Ptrial= 1; //Grado del polinomio (Ptrial = Ptest)
const unsigned int n_cycles=1;
FEMwDG<dim> Advection(Ptrial);
Advection.run(refinement_initial, typecase ,n_cycles, refine_type);
}
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;
}
return 0;
}
