Hello Wolfgang,
Yes, I am still getting the same error even using cell->set_future_fe_index(
0) before executing the refinement.
I have attached the error and all files needed to reproduce it, but
anything that concerns us is in main.cc.
the code is a mix between step-26 and step-27 to add elements during the
computation.
To do so, I set some elements with FE_Nothing and the beggining and then,
inside the method refine_mesh, I activate them at a given time_step using
set_future_fe_index as suggested by Bruno before executing
triangulation.execute_coarsening_and_refinement()
See:
if (timestep_number == 25) activate_FEs();
triangulation.execute_coarsening_and_refinement();
setup_system();
solution_trans.interpolate(previous_solution, solution);
when performing the interpolation with
solution_trans.interpolate(previous_solution,
solution); it gives me the error.
Thank you for the time,
El jueves, 18 de enero de 2024 a las 18:15:08 UTC, Wolfgang Bangerth
escribió:
> On 1/18/24 03:28, luis miguel reig buades wrote:
> >
> > However, it does not seem to solve the interpolation problem, I am
> assuming
> > the error happens because I did not assign an initial solution to the
> > activated elements.
> > Is there any way to assign the activated elements a solution as I am
> looping
> > through the cells and activating them? I have not been able to find this
>
> Are you saying that you continue to get the same error, even with the
> switched
> order of operations? Can you show us the code you use, and the error you
> get?
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/af9d8452-0e58-4e3f-ad22-8f36b94224d8n%40googlegroups.com.
// # physics constants
double global_PI = 3.1415927;
// Heat acitvation area
double square_length = 0.2;
// # Laser
double global_Pow_laser = 0.4; // power [W]
double global_spotsize_at_e_2 = 20e-6; // laser spot size at e^(-2) [m]
double global_c_laser = global_spotsize_at_e_2 / 4.0; // C parameter in Gaussian func, [m]
double global_c_hwhm = global_c_laser * 2.35482 / 2; // HWHM, [m]
// double global_V_scan_x = 10e-3; // scan speed, [m/s]
double global_V_scan_x = 4; // scan speed, [m/s]
double global_V_scan_y = 0; // scan speed, [m/s]
double global_init_position_x0 = 0; // initial spot center position
double global_init_position_y0 = 1;
// # material
// thin film
double global_rho_Tio2 = 4200; // mass density, [kg/m^3]
double global_C_Tio2 = 690; // heat capacity, [J/kg/K]
double global_k_Tio2 = 4.8; // thermal conductivity, [W/m/K]
// substrate
double global_rho_glass = 2200;
double global_C_glass = 700;
double global_k_glass = 1.8;
double global_film_thickness = 400e-9; // film thickness, [m]
// # simulation time
double global_simulation_time_step = 1e-5; // 10 [us]
double global_simulation_end_time = 100e-6 / global_V_scan_x; // 100 [um] / scan speed
// # about the MESH
#define BOUNDARY_NUM 11
/* ---------------------------------------------------------------------
*
* Copyright (C) 2013 - 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level output_directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, Texas A&M University, 2013
*/
// The program starts with the usual include files, all of which you should
// have seen before by now:
#include <deal.II/base/utilities.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/matrix_tools.h>
#include <fstream>
#include <iostream>
#include <filesystem>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/refinement.h>
// Then the usual placing of all content of this program into a namespace and
// the importation of the deal.II namespace into the one we will work in:
namespace Step26
{
using namespace dealii;
// Including the header file for the right hand side and boundary values
#ifndef GLOBAL_PARA
#define GLOBAL_PARA
#include "./globalPara.h"
// #include "./boundaryInit.h"
#include "./rightHandSide.h"
#endif
template <int dim>
class HeatEquation
{
public:
HeatEquation();
void run();
private:
void setup_system();
void solve_time_step();
void output_results() const;
void refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level);
// methods developed by me:
// void set_active_FEs();
void activate_FEs();
void deactivate_FEs();
void Create_Initial_Triangulation();
Triangulation<dim> triangulation;
const unsigned int initial_global_refinement;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> old_solution;
Vector<double> system_rhs;
double time;
double time_step;
unsigned int timestep_number;
const double theta;
std::string output_directory;
// attributes developed by me:
hp::FECollection<dim> fe_collection;
hp::QCollection<dim> quadrature_collection;
hp::QCollection<dim - 1> face_quadrature_collection;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim> & /*p*/,
const unsigned int component) const
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
template <int dim>
HeatEquation<dim>::HeatEquation()
: initial_global_refinement(2)
, fe(1)
, dof_handler(triangulation)
, time_step(1. / 500)
, theta(0.5)
, output_directory("output")
{
fe_collection.push_back(FE_Q<dim>(1));
fe_collection.push_back(FE_Nothing<dim>());
quadrature_collection.push_back(QGauss<dim>(2));
quadrature_collection.push_back(QGauss<dim>(2));
face_quadrature_collection.push_back(QGauss<dim - 1>(2));
face_quadrature_collection.push_back(QGauss<dim - 1>(2));
}
// @sect4{<code>HeatEquation::setup_system</code>}
//
// The next function is the one that sets up the DoFHandler object,
// computes the constraints, and sets the linear algebra objects
// to their correct sizes. We also compute the mass and Laplace
// matrix here by simply calling two functions in the library.
//
// Note that we do not take the hanging node constraints into account when
// assembling the matrices (both functions have an AffineConstraints argument
// that defaults to an empty object). This is because we are going to
// condense the constraints in run() after combining the matrices for the
// current time-step.
template <int dim>
void HeatEquation<dim>::setup_system()
{
dof_handler.distribute_dofs(fe_collection);
std::cout << std::endl
<< "===========================================" << std::endl
<< "Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl
<< std::endl;
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ true);
sparsity_pattern.copy_from(dsp);
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
system_matrix.reinit(sparsity_pattern);
MatrixCreator::create_mass_matrix(dof_handler,
quadrature_collection,
mass_matrix);
MatrixCreator::create_laplace_matrix(dof_handler,
quadrature_collection,
laplace_matrix);
solution.reinit(dof_handler.n_dofs());
old_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
// Method to specify the active and inactive FEs from the triangulation, which change every time powder is added.
template <int dim>
void HeatEquation<dim>::deactivate_FEs()
{
const Point<2> point(0, 0);
for (auto &cell: dof_handler.active_cell_iterators())
{
for (const auto v : cell->vertex_indices()){
const double distance_from_point =
point.distance(cell->vertex(v));
if (distance_from_point < 0.5){
cell->set_active_fe_index(1); // index 1 is for FE_Nothing elements, which are elements with 0 degrees of freedom, since it is the second element of the fe_collection array.
break;
}else{
cell->set_active_fe_index(0);
}
}
}
dof_handler.distribute_dofs(fe_collection);
}
// Method to activate all cells
template <int dim>
void HeatEquation<dim>::activate_FEs()
{
for (auto &cell: dof_handler.active_cell_iterators())
{
cell->set_future_fe_index(0);
}
dof_handler.distribute_dofs(fe_collection);
}
// @sect4{<code>HeatEquation::solve_time_step</code>}
//
// The next function is the one that solves the actual linear system
// for a single time step. There is nothing surprising here:
template <int dim>
void HeatEquation<dim>::solve_time_step()
{
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
SolverCG<Vector<double>> cg(solver_control);
PreconditionSSOR<SparseMatrix<double>> preconditioner;
preconditioner.initialize(system_matrix, 1.0);
cg.solve(system_matrix, solution, system_rhs, preconditioner);
constraints.distribute(solution);
std::cout << " " << solver_control.last_step() << " CG iterations."
<< std::endl;
}
// @sect4{<code>HeatEquation::output_results</code>}
//
// Neither is there anything new in generating graphical output other than the
// fact that we tell the DataOut object what the current time and time step
// number is, so that this can be written into the output file:
template <int dim>
void HeatEquation<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
//Output de degrees of freedom in every element:
Vector<float> fe_degrees(triangulation.n_active_cells());
for (const auto &cell : dof_handler.active_cell_iterators())
fe_degrees(cell->active_cell_index()) =
fe_collection[cell->active_fe_index()].degree;
data_out.add_data_vector(fe_degrees, "fe_degree");
data_out.add_data_vector(solution, "T");
data_out.build_patches();
data_out.set_flags(DataOutBase::VtkFlags(time, timestep_number));
// create output output_directory if it does not exist:
if (!std::filesystem::exists(output_directory)) (std::filesystem::create_directory(output_directory));
const std::string filename =
output_directory + "/solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk";
std::ofstream output(filename);
data_out.write_vtk(output);
}
// @sect4{<code>HeatEquation::refine_mesh</code>}
//
// This function is the interesting part of the program. It takes care of
// the adaptive mesh refinement. The three tasks
// this function performs is to first find out which cells to
// refine/coarsen, then to actually do the refinement and eventually
// transfer the solution vectors between the two different grids. The first
// task is simply achieved by using the well-established Kelly error
// estimator on the solution. The second task is to actually do the
// remeshing. That involves only basic functions as well, such as the
// <code>refine_and_coarsen_fixed_fraction</code> that refines those cells
// with the largest estimated error that together make up 60 per cent of the
// error, and coarsens those cells with the smallest error that make up for
// a combined 40 per cent of the error. Note that for problems such as the
// current one where the areas where something is going on are shifting
// around, we want to aggressively coarsen so that we can move cells
// around to where it is necessary.
//
// As already discussed in the introduction, too small a mesh leads to
// too small a time step, whereas too large a mesh leads to too little
// resolution. Consequently, after the first two steps, we have two
// loops that limit refinement and coarsening to an allowable range of
// cells:
template <int dim>
void HeatEquation<dim>::refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level)
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
face_quadrature_collection,
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
estimated_error_per_cell,
0.6,
0.4);
if (triangulation.n_levels() > max_grid_level)
for (const auto &cell :
triangulation.active_cell_iterators_on_level(max_grid_level))
cell->clear_refine_flag();
for (const auto &cell :
triangulation.active_cell_iterators_on_level(min_grid_level))
cell->clear_coarsen_flag();
// These two loops above are slightly different but this is easily
// explained. In the first loop, instead of calling
// <code>triangulation.end()</code> we may as well have called
// <code>triangulation.end_active(max_grid_level)</code>. The two
// calls should yield the same iterator since iterators are sorted
// by level and there should not be any cells on levels higher than
// on level <code>max_grid_level</code>. In fact, this very piece
// of code makes sure that this is the case.
// As part of mesh refinement we need to transfer the solution vectors
// from the old mesh to the new one. To this end we use the
// SolutionTransfer class and we have to prepare the solution vectors that
// should be transferred to the new grid (we will lose the old grid once
// we have done the refinement so the transfer has to happen concurrently
// with refinement). At the point where we call this function, we will
// have just computed the solution, so we no longer need the old_solution
// variable (it will be overwritten by the solution just after the mesh
// may have been refined, i.e., at the end of the time step; see below).
// In other words, we only need the one solution vector, and we copy it
// to a temporary object where it is safe from being reset when we further
// down below call <code>setup_system()</code>.
//
// Consequently, we initialize a SolutionTransfer object by attaching
// it to the old DoF handler. We then prepare the triangulation and the
// data vector for refinement (in this order).
SolutionTransfer<dim> solution_trans(dof_handler);
Vector<double> previous_solution;
previous_solution = solution;
triangulation.prepare_coarsening_and_refinement();
solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
// Now everything is ready, so do the refinement and recreate the DoF
// structure on the new grid, and finally initialize the matrix structures
// and the new vectors in the <code>setup_system</code> function. Next, we
// actually perform the interpolation of the solution from old to new
// grid. The final step is to apply the hanging node constraints to the
// solution vector, i.e., to make sure that the values of degrees of
// freedom located on hanging nodes are so that the solution is
// continuous. This is necessary since SolutionTransfer only operates on
// cells locally, without regard to the neighborhood.
if (timestep_number == 25) activate_FEs();
triangulation.execute_coarsening_and_refinement();
setup_system();
solution_trans.interpolate(previous_solution, solution);
constraints.distribute(solution);
}
template <int dim>
void HeatEquation<dim>::Create_Initial_Triangulation()
{
Point<dim> p1;
Point<dim> p2;
std::cout << "dim = " << dim << std::endl;
if (dim == 2){
p1 = {0, 0};
p2 = {2, 1};
} else if (dim == 3){
p1 = {0, 0, 0};
p2 = {2, 1, 1};
}
GridGenerator::hyper_rectangle(triangulation, p1 , p2);
triangulation.refine_global(initial_global_refinement);
}
template <int dim>
void HeatEquation<dim>::run()
{
const unsigned int n_adaptive_pre_refinement_steps = 4;
Create_Initial_Triangulation();
// set_active_FEs();
deactivate_FEs();
setup_system();
unsigned int pre_refinement_step = 0;
Vector<double> tmp;
Vector<double> forcing_terms;
start_time_iteration:
time = 0.0;
timestep_number = 0;
tmp.reinit(solution.size());
forcing_terms.reinit(solution.size());
VectorTools::interpolate(dof_handler,
Functions::ZeroFunction<dim>(),
old_solution);
solution = old_solution;
output_results();
// Then we start the main loop until the computed time exceeds our
// end time of 0.5. The first task is to build the right hand
// side of the linear system we need to solve in each time step.
// Recall that it contains the term $MU^{n-1}-(1-\theta)k_n AU^{n-1}$.
// We put these terms into the variable system_rhs, with the
// help of a temporary vector:
while (time <= 0.5)
{
time += time_step;
++timestep_number;
std::cout << "Time step " << timestep_number << " at t=" << time
<< std::endl;
mass_matrix.vmult(system_rhs, old_solution);
laplace_matrix.vmult(tmp, old_solution);
system_rhs.add(-(1 - theta) * time_step, tmp);
// The second piece is to compute the contributions of the source
// terms. This corresponds to the term $k_n
// \left[ (1-\theta)F^{n-1} + \theta F^n \right]$. The following
// code calls VectorTools::create_right_hand_side to compute the
// vectors $F$, where we set the time of the right hand side
// (source) function before we evaluate it. The result of this
// all ends up in the forcing_terms variable:
RightHandSide<dim> rhs_function;
rhs_function.set_time(time);
VectorTools::create_right_hand_side(dof_handler,
quadrature_collection,
rhs_function,
tmp);
forcing_terms = tmp;
forcing_terms *= time_step * theta;
rhs_function.set_time(time - time_step);
VectorTools::create_right_hand_side(dof_handler,
quadrature_collection,
rhs_function,
tmp);
forcing_terms.add(time_step * (1 - theta), tmp);
// Next, we add the forcing terms to the ones that
// come from the time stepping, and also build the matrix
// $M+k_n\theta A$ that we have to invert in each time step.
// The final piece of these operations is to eliminate
// hanging node constrained degrees of freedom from the
// linear system:
system_rhs += forcing_terms;
system_matrix.copy_from(mass_matrix);
system_matrix.add(theta * time_step, laplace_matrix);
constraints.condense(system_matrix, system_rhs);
// There is one more operation we need to do before we
// can solve it: boundary values. To this end, we create
// a boundary value object, set the proper time to the one
// of the current time step, and evaluate it as we have
// done many times before. The result is used to also
// set the correct boundary values in the linear system:
{
BoundaryValues<dim> boundary_values_function;
boundary_values_function.set_time(time);
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
boundary_values_function,
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
// With this out of the way, all we have to do is solve the
// system, generate graphical data, and...
solve_time_step();
output_results();
// ...take care of mesh refinement. Here, what we want to do is
// (i) refine the requested number of times at the very beginning
// of the solution procedure, after which we jump to the top to
// restart the time iteration, (ii) refine every fifth time
// step after that.
//
// The time loop and, indeed, the main part of the program ends
// with starting into the next time step by setting old_solution
// to the solution we have just computed.
if ((timestep_number == 1) &&
(pre_refinement_step < n_adaptive_pre_refinement_steps))
{
refine_mesh(initial_global_refinement,
initial_global_refinement +
n_adaptive_pre_refinement_steps);
++pre_refinement_step;
tmp.reinit(solution.size());
forcing_terms.reinit(solution.size());
std::cout << std::endl;
goto start_time_iteration;
}
else if ((timestep_number > 0) && (timestep_number % 5 == 0))
{
refine_mesh(initial_global_refinement,
initial_global_refinement +
n_adaptive_pre_refinement_steps);
tmp.reinit(solution.size());
forcing_terms.reinit(solution.size());
}
old_solution = solution;
}
}
} // namespace Step26
int main()
{
try
{
using namespace Step26;
HeatEquation<2> heat_equation_solver;
heat_equation_solver.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
--------------------------------------------------------
An error occurred in line <155> of file
</build/deal.ii-o0wOgt/deal.ii-9.5.1/source/dofs/dof_accessor_set.cc> in
function
void dealii::internal::set_dof_values(const dealii::DoFCellAccessor<dim,
spacedim, level_dof_access>&, const dealii::Vector<number>&, OutputVector&,
bool) [with int dim = 2; int spacedim = 2; bool lda = false; OutputVector =
dealii::Vector<double>; number = double]
The violated condition was:
local_values_old[i] == number() || get_abs(local_values_old[i] -
local_values[i]) <= get_abs(local_values_old[i] + local_values[i]) * 100000. *
std::numeric_limits<typename numbers::NumberTraits<
number>::real_type>::epsilon()
Additional information:
Called set_dof_values_by_interpolation(), but the element to be set,
value 0, does not match with the non-zero value 5.270338238394919e-05
already set before.
Stacktrace:
-----------
#0 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: void
dealii::internal::set_dof_values<2, 2, false, dealii::Vector<double>,
double>(dealii::DoFCellAccessor<2, 2, false> const&, dealii::Vector<double>
const&, dealii::Vector<double>&, bool)
#1 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: dealii::DoFCellAccessor<2, 2,
false>::set_dof_values_by_interpolation<dealii::Vector<double>,
double>(dealii::Vector<double> const&, dealii::Vector<double>&, unsigned short,
bool) const::{lambda(dealii::DoFCellAccessor<2, 2, false> const&,
dealii::Vector<double> const&,
dealii::Vector<double>&)#1}::operator()(dealii::DoFCellAccessor<2, 2, false>
const&, dealii::Vector<double> const&, dealii::Vector<double>&) const
#2 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: std::_Function_handler<void
(dealii::DoFCellAccessor<2, 2, false> const&, dealii::Vector<double> const&,
dealii::Vector<double>&), dealii::DoFCellAccessor<2, 2,
false>::set_dof_values_by_interpolation<dealii::Vector<double>,
double>(dealii::Vector<double> const&, dealii::Vector<double>&, unsigned short,
bool) const::{lambda(dealii::DoFCellAccessor<2, 2, false> const&,
dealii::Vector<double> const&,
dealii::Vector<double>&)#1}>::_M_invoke(std::_Any_data const&,
dealii::DoFCellAccessor<2, 2, false> const&, dealii::Vector<double> const&,
dealii::Vector<double>&)
#3 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: std::function<void
(dealii::DoFCellAccessor<2, 2, false> const&, dealii::Vector<double> const&,
dealii::Vector<double>&)>::operator()(dealii::DoFCellAccessor<2, 2, false>
const&, dealii::Vector<double> const&, dealii::Vector<double>&) const
#4 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: void
dealii::internal::process_by_interpolation<2, 2, false, dealii::Vector<double>,
double>(dealii::DoFCellAccessor<2, 2, false> const&, dealii::Vector<double>
const&, dealii::Vector<double>&, unsigned short, std::function<void
(dealii::DoFCellAccessor<2, 2, false> const&, dealii::Vector<double> const&,
dealii::Vector<double>&)> const&)
#5 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: void
dealii::DoFCellAccessor<2, 2,
false>::set_dof_values_by_interpolation<dealii::Vector<double>,
double>(dealii::Vector<double> const&, dealii::Vector<double>&, unsigned short,
bool) const
#6 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: dealii::SolutionTransfer<2,
dealii::Vector<double>, 2>::interpolate(std::vector<dealii::Vector<double>,
std::allocator<dealii::Vector<double> > > const&,
std::vector<dealii::Vector<double>, std::allocator<dealii::Vector<double> > >&)
const
#7 /lib/x86_64-linux-gnu/libdeal.ii.g.so.9.5.1: dealii::SolutionTransfer<2,
dealii::Vector<double>, 2>::interpolate(dealii::Vector<double> const&,
dealii::Vector<double>&) const
#8 ./main: Step26::HeatEquation<2>::refine_mesh(unsigned int, unsigned int)
#9 ./main: Step26::HeatEquation<2>::run()
#10 ./main: main
--------------------------------------------------------
make[3]: *** [CMakeFiles/run.dir/build.make:58: CMakeFiles/run] Aborted (core
dumped)
make[2]: *** [CMakeFiles/Makefile2:218: CMakeFiles/run.dir/all] Error 2
make[1]: *** [CMakeFiles/Makefile2:225: CMakeFiles/run.dir/rule] Error 2
make: *** [Makefile:183: run] Error 2//#----------------------------------------------------------
//#
//# This file defines the boundary and initial conditions
//#
//#----------------------------------------------------------
#ifndef GLOBAL_PARA
#define GLOBAL_PARA
#include "./globalPara.h"
#endif
//# Declaration
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide() : Function<dim>() {}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
//# Implementation
template <int dim>
double RightHandSide<dim>::value(const Point<dim> & p,
const unsigned int component) const
{
(void)component;
AssertIndexRange(component, 1);
Assert(dim == 2, ExcNotImplemented());
const double time = this->get_time();
double position_square_x = global_init_position_x0 + global_V_scan_x * time;
double position_square_y = global_init_position_y0 + global_V_scan_y * time;
if ( (p[0] > position_square_x - square_length) && (p[0] < position_square_x + square_length)
&&
(p[1] > position_square_y - square_length) && (p[1] < position_square_y + square_length)
){
return 2;
}
else{
return 0;
}
}
