Hello Daniel,
Please find the attached .cc file along with the .txt file. The attached
program runs perfectly
and gives a rectangular box. Now when I want to create the parallelepiped
by uncommenting out
the commented portions in 'grid_y_transform' and in 'system_setup'
functions (of course the
corresponding lines in the present file have to be commented then), I get
the following error (in red).
/home/abasak/Desktop/Anup/dealiicodes/myprogs/adaptime_subsequent_eta/half_sample/to_Daniel/PhaseField.cc:237:84:
error: no matching function for call to ‘dealii::Tensor<1, 3,
double>::Tensor(double, double, double)’
Tensor<1, dim>(0., 0.,
width*tan_theta*scale ));
Thanks a lot.
Anup.
On Tuesday, November 29, 2016 at 10:32:36 AM UTC-6, Daniel Arndt wrote:
>
> Anup,
>
> If you provide something compilable, I can have a look. ;-)
>
> Best,
> Daniel
>
> Am Dienstag, 29. November 2016 17:22:39 UTC+1 schrieb Anup Basak:
>>
>> Hello Daniel,
>>
>> I am attaching a slightly modified version for better readability. In
>> actual code
>> I read the values from a .prm file.
>>
>> Thanks,
>> Anup.
>>
>> On Tuesday, November 29, 2016 at 9:52:40 AM UTC-6, Anup Basak wrote:
>>>
>>> Hello Daniel,
>>>
>>> Thanks for your reply and help. I am attaching a part of the code
>>> herewith.
>>>
>>> Thanks and regards,
>>> Anup.
>>>
>>> On Tuesday, November 29, 2016 at 9:41:19 AM UTC-6, Daniel Arndt wrote:
>>>>
>>>> Anup,
>>>>
>>>> I am trying to implement the periodic boundary condition for a
>>>>> nonlinear parabolic equation (scalar).
>>>>> I have considered a *parallelepiped *as shown in the attached pdf. I
>>>>> want to impose periodic
>>>>> BC at x=0 and x=w surfaces. I have used an offset vector as appearing
>>>>> in the piece of code
>>>>> below, which I formed based on my understanding from the documentation
>>>>> (it is tangential to the
>>>>> surface x = w). Could you please tell
>>>>> me if it is the correct expression. I think the problem is in this
>>>>> vector only, because when I reduce the
>>>>> angle theta to zero, it works. Boundary x =0 and x=w are numbered as 1
>>>>> and 11, respectively.
>>>>> I also tried with the offset vector : Tensor<1, dim>(0., 0., -
>>>>> w*tan(theta))
>>>>> Any suggestion and help would be appreciated.
>>>>>
>>>>> constraints_eta1.clear();
>>>>> std::vector<GridTools::PeriodicFacePair<typename
>>>>> DoFHandler<dim>::cell_iterator> > periodicity_vector1;
>>>>> GridTools::collect_periodic_faces(dof_handler_eta, 1, 11, 0,
>>>>> periodicity_vector1,
>>>>> Tensor<1, dim>(0., 0.,
>>>>> w*tan(theta)));
>>>>>
>>>>> DoFTools::make_periodicity_constraints<DoFHandler<dim>
>>>>> >(periodicity_vector1, constraints_eta1);
>>>>> constraints_eta1.close();
>>>>>
>>>> That looks and sounds correct if the vertices of the bounding
>>>> parallelepiped have the coordinates
>>>> (0|0|0), (w|0|tan(theta)*w), (0|h|0), (w|h|tan(theta)*w), (0|0|L),
>>>> (w|0|L+tan(theta)*w), (0|h|L), (w|h|L+tan(theta)*w)
>>>>
>>>> Can you confirm that the coarse cells have the correct boundary
>>>> indicators? Would you have a minimal code?
>>>>
>>>> Best,
>>>> Daniel
>>>>
>>>
--
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].
For more options, visit https://groups.google.com/d/optout.
##
# CMake script for the step-44 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "PhaseField")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.4 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
#include <deal.II/base/function.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/std_cxx11/shared_ptr.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/fe_dgp_monomial.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/compressed_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition_selector.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/matrix_out.h>
#include <iostream>
#include <fstream>
//////////////////////////////////////////////////////////////////
namespace PhaseField
{
using namespace dealii;
///////////////////////////////////
template <int dim>
class Solid
{
public:
Solid();
virtual
~Solid();
void
run();
void
make_grid(); // function creates geometry and mesh
void
system_setup(); // sets of all the vectors and matrices for Ginzburg-Landau eqn and the mechanics problem
void
output_results() const;
Triangulation<dim> triangulation;
const unsigned int degree_eta; // degree of polynomial for eta
FE_Q<dim> fe_eta; // fe object for eta
DoFHandler<dim> dof_handler_eta; //another dof_handler for eta
ConstraintMatrix constraints_eta1;
SparsityPattern sparsity_pattern_eta1;
SparseMatrix<double> mass_matrix1; // mass matrix out of Ginxburg-Landau eqn
SparseMatrix<double> laplace_matrix1; // Laplace matrix out of Ginxburg-Landau eqn
SparseMatrix<double> system_matrix1; // Laplace matrix out of Ginxburg-Landau eqn
SparseMatrix<double> nl_matrix1; //
Vector<double> solution_eta1; // solution vector for eta
Vector<double> solution_eta2; // solution vector for eta
};
/////////////////////////////////////////////////////////
// defines the initial condition for the order parameter
template <int dim>
class Initial_eta1 : public Function<dim>
{
public:
Initial_eta1 () : Function<dim>() {}
virtual double value(const Point<dim> &p,
const unsigned int /*component = 0*/) const;
};
template <int dim>
double Initial_eta1<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
if (p[0] <= 5.25e-9 && p[2] <= 30.0e-9+0.54138*p[0])
return 1.0;
else if ((p[0] > 5.25e-9 && p[2] <= 30.0e-9+0.54138*p[0]) &&
(p[0] <= 9.75e-9 && p[2] <= 30.0e-9+0.54138*p[0]))
return 0.0;
else if (p[0] > 9.75e-9 && p[2] <= 30.0e-9+0.54138*p[0])
return 1.0;
else return 0.0;
}
/////////////////////////////////////////////////////////
// defines the initial condition for the order parameter
template <int dim>
class Initial_eta2 : public Function<dim>
{
public:
Initial_eta2 () : Function<dim>() {}
virtual double value(const Point<dim> &p,
const unsigned int /*component = 0*/) const;
};
template <int dim>
double Initial_eta2<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
if (p[0] <= 5.25e-9 && p[2] <= 30.0e-9+0.54138*p[0])
return 0.0;
else if ((p[0] > 5.25e-9 && p[2] <= 30.0e-9+0.54138*p[0]) &&
(p[0] <= 9.75e-9 && p[2] <= 30.0e-9+0.54138*p[0]))
return 1.0;
else if (p[0] > 9.75e-9 && p[2] <= 30.0e-9+0.54138*p[0])
return 0.0;
else return 0.0;
}
///////////////////////////////////////////////////
template <int dim>
Solid<dim>::Solid()
:
triangulation(Triangulation<dim>::maximum_smoothing),
degree_eta(1),
fe_eta (1),
dof_handler_eta (triangulation)
{
}
template <int dim>
Solid<dim>::~Solid()
{
dof_handler_eta.clear();
}
////////////////// current /////////////////////
template <int dim>
void Solid<dim>::run()
{
make_grid();
system_setup();
VectorTools::interpolate(dof_handler_eta, Initial_eta1<dim>(), solution_eta1);
VectorTools::interpolate(dof_handler_eta, Initial_eta2<dim>(), solution_eta2);
output_results();
}
///////////////////////////////////////////////////////
template <int dim>
Point<dim> grid_y_transform (const Point<dim> &pt_in)
{
const double &x = pt_in[0];
const double &z = pt_in[2];
/*
const double z_upper = 45.0 + (5.413861341550015/10.0)*x; // Line defining upper edge of beam
const double z_lower = 0.0 + (5.413861341550015/10.0)*x; // Line defining lower edge of beam
const double theta = z/45.0; // Fraction of height along left side of beam
const double y_transform = (1-theta)*z_lower + theta*z_upper; // Final transformation
*/
const double z_upper = 45.0 ; // Line defining upper edge of beam
const double z_lower = 0.0; // Line defining lower edge of beam
const double theta = z/45.0; // Fraction of height along left side of beam
const double y_transform = (1-theta)*z_lower + theta*z_upper; // Final transformation
Point<dim> pt_out = pt_in;
pt_out[2] = y_transform;
return pt_out;
}
template <int dim>
void Solid<dim>::make_grid()
{
std::vector< unsigned int > repetitions(dim, 135);
repetitions[dim-3] = 45;
if (dim == 3)
repetitions[dim-2] = 1;
const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0));
const Point<dim> top_right = (dim == 3 ? Point<dim>(15, 2, 45) :
Point<dim>(15, 45));
GridGenerator::subdivided_hyper_rectangle(triangulation,
repetitions,
bottom_left,
top_right);
const double tol_boundary = 1e-6;
typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(), endc = triangulation.end();
for (; cell != endc; ++cell)
for (unsigned int face = 0;
face < GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary() == true)
{
if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
cell->face(face)->set_boundary_id(1); // -X faces
else if (std::abs(cell->face(face)->center()[0] - 15) < tol_boundary)
cell->face(face)->set_boundary_id(11); // +X faces
}
GridTools::transform(&grid_y_transform<dim>, triangulation);
GridTools::scale(1e-9, triangulation);
}
///////////////////////////////////////////////////
template <int dim>
void Solid<dim>::system_setup()
{
dof_handler_eta.distribute_dofs (fe_eta);
const unsigned int n_dofs_eta = dof_handler_eta.n_dofs();
double width = 15;
double scale = 1.0e-9;
double tan_theta = 0.5413861341550015;
constraints_eta1.clear();
std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> > periodicity_vector1;
GridTools::collect_periodic_faces(dof_handler_eta, 1, 11, 0, periodicity_vector1,
Tensor<1, dim>());
/* GridTools::collect_periodic_faces(dof_handler_eta, 1, 11, 0, periodicity_vector1,
Tensor<1, dim>(0., 0., width*tan_theta*scale ));
*/
DoFTools::make_periodicity_constraints<DoFHandler<dim> >(periodicity_vector1, constraints_eta1);
constraints_eta1.close();
{
mass_matrix1.clear ();
laplace_matrix1.clear ();
system_matrix1.clear ();
nl_matrix1.clear ();
DynamicSparsityPattern dsp1 (n_dofs_eta, n_dofs_eta);
DoFTools::make_sparsity_pattern (dof_handler_eta, dsp1, constraints_eta1, true);
sparsity_pattern_eta1.copy_from(dsp1);
mass_matrix1.reinit (sparsity_pattern_eta1);
laplace_matrix1.reinit (sparsity_pattern_eta1);
system_matrix1.reinit (sparsity_pattern_eta1);
nl_matrix1.reinit (sparsity_pattern_eta1);
MatrixCreator::create_mass_matrix(dof_handler_eta,
QGauss<dim>(fe_eta.degree+1), //creates the mass matrix
mass_matrix1);
MatrixCreator::create_laplace_matrix(dof_handler_eta,
QGauss<dim>(fe_eta.degree+1), //creates the Laplace matrix
laplace_matrix1);
solution_eta1.reinit(dof_handler_eta.n_dofs());
solution_eta2.reinit(dof_handler_eta.n_dofs());
}
}
///////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler_eta);
data_out.add_data_vector(solution_eta1, "eta1");
data_out.build_patches();
const std::string filename = "solution-"
+ Utilities::int_to_string(10, 3) +
".vtk";
std::ofstream output(filename.c_str());
data_out.write_vtk(output);
}
}
/////////////////////////////////////////////////////////
int main ()
{
using namespace dealii;
using namespace PhaseField;
try
{
deallog.depth_console(0);
Solid<3> solid_3d;
solid_3d.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl << exc.what()
<< std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl << "Aborting!"
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
