Luca,
Thank you for getting back to me so quickly. Unfortunately I am still
stuck.
I sort of understand what you are saying here:
> it should be sufficient to copy the class you are interested in and make
> it derived from FlatManifold instead of Boundary
but I don't know exactly which class you are referring to or how to
implement that. My guess is something like:
static OpenCASCADE::NormalProjectionBoundary<2,3>
normal_projector(cad_surface, tolerance);
FlatManifold<2,3> normal_projector_as_FlatManifold = FlatManifold(
*normal_projector_copy*);
triangulation.set_manifold(1,normal_projector_as_FlatManifold);
where I am confused as to how to obtain *normal_projector_copy*. I'm also
uncertain if this is the class(?) you are referring to.
Additionally, you provided the helpful remark:
> By the way, you are hitting this problem because you try to use higher
> order mappings on co-dimension one surfaces, defined through iges files
I see that I am using the default degree=2 to define the order of my
mappings on the codimension one surface. I ran the laplace_on_surface code
using degree = 1 and this moved my error to another location in the code
base! This is a positive sign, but doesn't get me all the way home. I've
attached a modified laplace_on_surface.cc file which has the following two
features:
1. an option in make_grid_and_dofs() which specifies where to get the
mesh and manifold description from. The string manifold_method specifies
"cad" or "GridGenerator", and I have implemented a similar codimension one
sphere from the GridGenerator utility.
2. I have supplied the degree argument in the instantiation of the
LaplaceProblem<spacedim> class in the main() function, this currently is
set to 1 (rather than the default, degree=2).
When I run the new function on the sphere generated using GridGenerator
life is good (for degree 1 and 2). When I run the new function with the
sphere generated by manifold_method="cad" and degree 2 I get the same old
problem (as expected) - when I run with degree 1, the old problem appears
to go away but the solver does not converge.
> ----------------------------------------------------
> Exception on processing:
> --------------------------------------------------------
> An error occurred in line <606> of file
> </usr/local/dealii/include/deal.II/lac/solver_cg.h> in function
> void dealii::SolverCG<VectorType>::solve(const MatrixType&,
> VectorType&, const VectorType&, const PreconditionerType&) [with MatrixType
> = dealii::SparseMatrix<double>; PreconditionerType =
> dealii::PreconditionSSOR<>; VectorType = dealii::Vector<double>]
> The violated condition was:
> false
> The name and call sequence of the exception was:
> SolverControl::NoConvergence (it, res)
> Additional Information:
> Iterative method reported convergence failure in step 3074. The residual
> in the last step was 2.09928.
This is unexpected because the sphere generated from cad should be
basically the same as the sphere generated by GridGenerator. It seems to
me that I still cannot treat meshes and manifolds brought in through cad as
though they were the same as ones generated by GridGenerator.
Thank you for your time,
Tom
On Wed, Jun 29, 2016 at 6:39 AM, luca.heltai <[email protected]> wrote:
> Dear Tom,
>
> you’re not doing anything wrong. The Manifold/Boundary interfaces are
> undergoing a lot of restructuring, and you hit a problem of backward
> compatibility implementation…
>
> Your issue is in the following lines of code, that were inserted to
> maintain backward compatibility between Manifold and Boundary objects:
>
>
> https://github.com/dealii/dealii/pull/2418/files#diff-c743c51b4f2c95b842d938d7f90da390L3545
>
> The above pull request is work in progress, as it breaks a lot of backward
> compatibilities. I hope to have some time soon to work on it...
>
> In particular, you see that if an object is derived by a Boundary class,
> then a different code is executed w.r.t pure Manifold classes (Boundary is
> derived from Manifold, but knows of being a codimension one surface).
>
> Unfortunately the OpenCascade wrappers were derived from Boundary (they
> were introduced before the Manifold concept was formulated), but they do
> not implement all function calls
>
> get_intermediate_points_on_line/quad/hex
>
> which is what triggers your exception. In your case it should be
> sufficient to copy the class you are interested in and make it derived from
> FlatManifold instead of
> Boundary. Everything should then work out of the box.
>
> By the way, you are hitting this problem because you try to use higher
> order mappings on co-dimension one surfaces, defined through iges files.
> This should works, but its highly untested… you may hit some high order
> mapping bugs we are trying to address… I hope things will be resolved soon.
>
> Can you open an issue, so that I can refer back to it when things are
> resolved?
>
> Best,
>
> Luca.
>
>
> > On 29 Jun 2016, at 24:19, thomas stephens <[email protected]> wrote:
> >
> > I am trying to combine techniques from tutorial 54 and tutorial 38 -
> namely, I would like to solve a surface diffusion problem on a geometry
> that I bring in through the opencascade interface.
> >
> > I have created a geometry (using Trelis/Cubit) and have the mesh .ucd
> file, as well as the .iges file (attached as sphere_mesh.ucd and
> sphere.iges). I can build a GridIn<2,3> object from sphere_mesh.ucd and
> attach that to a Triangulation, and I can create a TopoDS_Shape object from
> sphere.iges, and I can attach a NormalProjectionBoundary<2,3>
> normal_projector to the triangulation - at least to the extent that I can
> then call triangulation.refine_global(n) and obtain a refinement of
> sphere_mesh.ucd that continues to yield better approximations to a sphere.
> (This all occurs in the make_grid_and_dofs() function in
> laplace_on_surface.cc)
> >
> > The problem I am having shows up in the assemble_system() function in
> laplace_on_surface.cc. The error I get is
> >
> >
> > An error occurred in line <53> of file
> </home/tds3/software/dealii/dealii-8.4.1/source/grid/tria_boundary.cc> in
> function
> > void dealii::Boundary<dim,
> spacedim>::get_intermediate_points_on_quad(const typename
> dealii::Triangulation<dim, spacedim>::quad_iterator&,
> std::vector<dealii::Point<spacedim> >&) const [with int dim = 2; int
> spacedim = 3; typename dealii::Triangulation<dim, spacedim>::quad_iterator
> = dealii::TriaIterator<dealii::CellAccessor<2, 3> >]
> > The violated condition was:
> > false
> > The name and call sequence of the exception was:
> > ExcPureFunctionCalled()
> > You (or a place in the library) are trying to call a function that is
> declared as a virtual function in a base class but that has not been
> overridden in your derived class.
> >
> > The entire stacktrace is:
> > #0 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::Boundary<2,
> 3>::get_intermediate_points_on_quad(dealii::TriaIterator<dealii::CellAccessor<2,
> 3> > const&, std::vector<dealii::Point<3, double>,
> std::allocator<dealii::Point<3, double> > >&) const
> > #1 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1:
> > #2 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1:
> dealii::MappingQGeneric<2,
> 3>::add_quad_support_points(dealii::TriaIterator<dealii::CellAccessor<2, 3>
> > const&, std::vector<dealii::Point<3, double>,
> std::allocator<dealii::Point<3, double> > >&) const
> > #3 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1:
> dealii::MappingQGeneric<2,
> 3>::compute_mapping_support_points(dealii::TriaIterator<dealii::CellAccessor<2,
> 3> > const&) const
> > #4 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1:
> dealii::MappingQGeneric<2,
> 3>::fill_fe_values(dealii::TriaIterator<dealii::CellAccessor<2, 3> >
> const&, dealii::CellSimilarity::Similarity, dealii::Quadrature<2> const&,
> dealii::Mapping<2, 3>::InternalDataBase const&,
> dealii::internal::FEValues::MappingRelatedData<2, 3>&) const
> > #5 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::MappingQ<2,
> 3>::fill_fe_values(dealii::TriaIterator<dealii::CellAccessor<2, 3> >
> const&, dealii::CellSimilarity::Similarity, dealii::Quadrature<2> const&,
> dealii::Mapping<2, 3>::InternalDataBase const&,
> dealii::internal::FEValues::MappingRelatedData<2, 3>&) const
> > #6 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: dealii::FEValues<2,
> 3>::do_reinit()
> > #7 /usr/local/dealii/lib/libdeal_II.g.so.8.4.1: void
> dealii::FEValues<2, 3>::reinit<dealii::DoFHandler,
> false>(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2,
> 3>, false> > const&)
> > #8 ./laplace_on_surface:
> laplace_on_surface::LaplaceBeltramiProblem<3>::assemble_system()
> > #9 ./laplace_on_surface:
> laplace_on_surface::LaplaceBeltramiProblem<3>::run()
> > #10 ./laplace_on_surface: main
> >
> > It seems to me that I am missing something in the way I bring in the
> geometry back in make_grid_and_dofs() because I can generate a spherical
> mesh using the GridGenerator objects (as is done in tutorial step 38, and
> everything works fine. My suspicion is that I am not setting manifold ids
> correctly in my current implementation above.
> >
> > Also, please suggest appropriate tags for this question - I thought this
> should get an opencascade tag, but that's not a 'suggested tag'.
> >
> >
> > Thank you,
> > Tom
> >
> > --
> > 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.
> > <CMakeLists.txt><laplace_on_surface.cc><sphere_mesh.ucd><sphere.iges>
>
> --
> 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 a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/XG6J8oX1XxQ/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> [email protected].
> For more options, visit https://groups.google.com/d/optout.
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
/*
* Copyright (C) 2010 - 2015 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Authors: Andrea Bonito, Sebastian Pauletti.
*/
// @sect3{Include files}
// If you've read through step-4 and step-7, you will recognize that we have
// used all of the following include files there already. Consequently, we
// will not explain their meaning here again.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.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/opencascade/boundary_lib.h>
#include <deal.II/opencascade/utilities.h>
#include <fstream>
#include <iostream>
namespace laplace_on_surface
{
using namespace dealii;
//// class LaplaceBeltramiProblem
template <int spacedim>
class LaplaceBeltramiProblem
{
public:
LaplaceBeltramiProblem (const unsigned degree = 2);
void run ();
private:
static const unsigned int dim = spacedim-1;
void make_grid_and_dofs ();
void assemble_system ();
void solve ();
void output_results () const;
void compute_error () const;
Triangulation<dim,spacedim> triangulation;
FE_Q<dim,spacedim> fe;
DoFHandler<dim,spacedim> dof_handler;
MappingQ<dim, spacedim> mapping;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
//// class Solution
template <int dim>
class Solution : public Function<dim>
{
public:
Solution () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual Tensor<1,dim> gradient (const Point<dim> &p,
const unsigned int component = 0) const;
};
// function definitions for Solution class
/*{{{*/
template <>
double Solution<3>::value (const Point<3> &p, const unsigned int) const {
return (std::sin(numbers::PI * p(0)) *
std::cos(numbers::PI * p(1))*exp(p(2)));
}
template <>
Tensor<1,3> Solution<3>::gradient (const Point<3> &p, const unsigned int) const {
using numbers::PI;
Tensor<1,3> return_value;
return_value[0] = PI *cos(PI * p(0))*cos(PI * p(1))*exp(p(2));
return_value[1] = -PI *sin(PI * p(0))*sin(PI * p(1))*exp(p(2));
return_value[2] = sin(PI * p(0))*cos(PI * p(1))*exp(p(2));
return return_value;
}
/*}}}*/
// class RightHandSide
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;
};
// function definitions for RightHandSide class
/*{{{*/
template <>
double RightHandSide<3>::value (const Point<3> &p, const unsigned int /*component*/) const {
/*{{{*/
using numbers::PI;
Tensor<2,3> hessian;
hessian[0][0] = -PI*PI*sin(PI*p(0))*cos(PI*p(1))*exp(p(2));
hessian[1][1] = -PI*PI*sin(PI*p(0))*cos(PI*p(1))*exp(p(2));
hessian[2][2] = sin(PI*p(0))*cos(PI*p(1))*exp(p(2));
hessian[0][1] = -PI*PI*cos(PI*p(0))*sin(PI*p(1))*exp(p(2));
hessian[1][0] = -PI*PI*cos(PI*p(0))*sin(PI*p(1))*exp(p(2));
hessian[0][2] = PI*cos(PI*p(0))*cos(PI*p(1))*exp(p(2));
hessian[2][0] = PI*cos(PI*p(0))*cos(PI*p(1))*exp(p(2));
hessian[1][2] = -PI*sin(PI*p(0))*sin(PI*p(1))*exp(p(2));
hessian[2][1] = -PI*sin(PI*p(0))*sin(PI*p(1))*exp(p(2));
Tensor<1,3> gradient;
gradient[0] = PI * cos(PI*p(0))*cos(PI*p(1))*exp(p(2));
gradient[1] = - PI * sin(PI*p(0))*sin(PI*p(1))*exp(p(2));
gradient[2] = sin(PI*p(0))*cos(PI*p(1))*exp(p(2));
Point<3> normal = p;
normal /= p.norm();
return (- trace(hessian)
+ 2 * (gradient * normal)
+ (hessian * normal) * normal);
/*}}}*/
}
/*}}}*/
// constructor for LaplaceBeltramiOperator
template <int spacedim>
LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem (const unsigned degree)
:
fe (degree),
dof_handler(triangulation),
mapping (degree)
{}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs () {
/*{{{*/
const std::string manifold_method = "cad";
//const std::string manifold_method = "GridGenerator";
if (manifold_method == "cad") {
cout << "using cad generated mesh" << endl;
const std::string cad_file_name = "sphere.iges";
TopoDS_Shape cad_surface = OpenCASCADE::read_IGES(cad_file_name, 1);
const double tolerance = OpenCASCADE::get_shape_tolerance(cad_surface) * 5;
std::vector<TopoDS_Compound> compounds;
std::vector<TopoDS_CompSolid> compsolids;
std::vector<TopoDS_Solid> solids;
std::vector<TopoDS_Shell> shells;
std::vector<TopoDS_Wire> wires;
OpenCASCADE::extract_compound_shapes(cad_surface,
compounds,
compsolids,
solids,
shells,
wires);
std::ifstream in;
std::string in_mesh_filename = "sphere_mesh.ucd";
in.open(in_mesh_filename.c_str());
GridIn<2,3> gi;
gi.attach_triangulation(triangulation);
gi.read (in);
Triangulation<2,3>::active_cell_iterator cell = triangulation.begin_active();
cell->set_all_manifold_ids(1);
Assert(wires.size() > 0,
ExcMessage("I could not find any wire in the CAD file you gave me. Bailing out."));
static OpenCASCADE::NormalProjectionBoundary<2,3> normal_projector(cad_surface, tolerance);
triangulation.set_manifold(1,normal_projector);
triangulation.refine_global(3);
}
else if (manifold_method == "GridGenerator") {
cout << "using GridGenerator generated mesh" << endl;
static SphericalManifold<dim,spacedim> surface_description;
GridGenerator::hyper_sphere(triangulation);
triangulation.set_all_manifold_ids(1);
triangulation.set_manifold (1, surface_description);
triangulation.refine_global(3);
}
else {
cout << "in make_grid_and_dofs(): manifold_method not set correctly" << endl;
exit(1);
}
// output results
const std::string out_filename = "sphere_surface.vtk";
std::ofstream logfile(out_filename.c_str());
GridOut grid_out;
grid_out.write_vtk(triangulation, logfile);
cout << "output file written" << endl;
//
std::cout << "Surface mesh has " << triangulation.n_active_cells()
<< " active cells."
<< std::endl;
dof_handler.distribute_dofs (fe);
std::cout << "Surface mesh has " << dof_handler.n_dofs()
<< " degrees of freedom."
<< std::endl;
DynamicSparsityPattern dsp (dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
sparsity_pattern.copy_from (dsp);
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
/*}}}*/
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::assemble_system () {
/*{{{*/
system_matrix = 0;
system_rhs = 0;
const QGauss<dim> quadrature_formula(2*fe.degree);
FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula,
update_values |
update_gradients |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<double> rhs_values(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
const RightHandSide<spacedim> rhs;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
cell!=endc; ++cell)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
rhs.value_list (fe_values.get_quadrature_points(), rhs_values);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_matrix(i,j) += fe_values.shape_grad(i,q_point) *
fe_values.shape_grad(j,q_point) *
fe_values.JxW(q_point);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_rhs(i) += fe_values.shape_value(i,q_point) *
rhs_values[q_point]*
fe_values.JxW(q_point);
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
//std::map<types::global_dof_index,double> boundary_values;
//VectorTools::interpolate_boundary_values (mapping,
// dof_handler,
// 0,
// Solution<spacedim>(),
// boundary_values);
//MatrixTools::apply_boundary_values (boundary_values,
// system_matrix,
// solution,
// system_rhs,false);
/*}}}*/
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::solve ()
{
/*{{{*/
SolverControl solver_control (solution.size(),
1e-7 * system_rhs.l2_norm());
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
/*}}}*/
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::output_results () const
{
/*{{{*/
DataOut<dim,DoFHandler<dim,spacedim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution,
"solution",
DataOut<dim,DoFHandler<dim,spacedim> >::type_dof_data);
data_out.build_patches (mapping,
mapping.get_degree());
std::string filename ("solution-");
filename += static_cast<char>('0'+spacedim);
filename += "d.vtk";
std::ofstream output (filename.c_str());
data_out.write_vtk (output);
/*}}}*/
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::compute_error () const
{
/*{{{*/
Vector<float> difference_per_cell (triangulation.n_active_cells());
VectorTools::integrate_difference (mapping, dof_handler, solution,
Solution<spacedim>(),
difference_per_cell,
QGauss<dim>(2*fe.degree+1),
VectorTools::H1_norm);
std::cout << "H1 error = "
<< difference_per_cell.l2_norm()
<< std::endl;
/*}}}*/
}
template <int spacedim>
void LaplaceBeltramiProblem<spacedim>::run () {
/*{{{*/
cout << "entered run ()" << endl;
make_grid_and_dofs();
cout << "made grid and dofs" << endl;
assemble_system ();
cout << "assembled system" << endl;
solve ();
cout << "solved" << endl;
output_results ();
compute_error ();
/*}}}*/
}
}
// @sect3{The main() function}
// The remainder of the program is taken up by the <code>main()</code>
// function. It follows exactly the general layout first introduced in step-6
// and used in all following tutorial programs:
int main ()
{
try
{
using namespace dealii;
using namespace laplace_on_surface;
const unsigned int spacedim = 3;
LaplaceBeltramiProblem<spacedim> laplace_beltrami(1); // default: degree=2
cout << "about to run laplace_beltrami..." << endl;
laplace_beltrami.run();
cout << "ran laplace_beltrami..." << endl;
}
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;
}