Commit d6540723 authored by simone's avatar simone
Browse files

taking in to account changes to dataMesh, where mesh is not protected anymore. Now it is private.

parent 9bd984b0
......@@ -190,12 +190,12 @@ ElasticStructureHandler( const GetPot& data_file,
BCHandler& BCh ):
DataElasticStructure<Mesh>( data_file ),
_refFE( refFE ),
_dof( this->_mesh, _refFE ),
_dof( this->mesh(), _refFE ),
_dim( _dof.numTotalDof() ),
_Qr( Qr ),
_bdQr( bdQr ),
_fe( _refFE, getGeoMap( this->_mesh ), _Qr ),
_feBd( _refFE.boundaryFE(), getGeoMap( this->_mesh ).boundaryMap(), _bdQr ),
_fe( _refFE, getGeoMap( this->mesh() ), _Qr ),
_feBd( _refFE.boundaryFE(), getGeoMap( this->mesh() ).boundaryMap(), _bdQr ),
_d( _dim ),
_dRhs( _dim ),
_w( _dim ),
......@@ -212,12 +212,12 @@ ElasticStructureHandler( const GetPot& data_file,
const QuadRule& bdQr):
DataElasticStructure<Mesh>( data_file ),
_refFE( refFE ),
_dof( this->_mesh, _refFE ),
_dof( this->mesh(), _refFE ),
_dim( _dof.numTotalDof() ),
_Qr( Qr ),
_bdQr( bdQr ),
_fe( _refFE, getGeoMap( this->_mesh ), _Qr ),
_feBd( _refFE.boundaryFE(), getGeoMap( this->_mesh ).boundaryMap(), _bdQr ),
_fe( _refFE, getGeoMap( this->mesh() ), _Qr ),
_feBd( _refFE.boundaryFE(), getGeoMap( this->mesh() ).boundaryMap(), _bdQr ),
_d( _dim ),
_dRhs( _dim ),
_w( _dim ),
......@@ -290,25 +290,25 @@ ElasticStructureHandler<Mesh>::postProcess()
}
namedef = "defor." + name + ".mesh";
wr_medit_ascii_scalar( "dep_x." + name + ".bb", _d.giveVec(), this->_mesh.numVertices() );
wr_medit_ascii_scalar( "dep_y." + name + ".bb", _d.giveVec() + _dim, this->_mesh.numVertices() );
wr_medit_ascii_scalar( "dep_z." + name + ".bb", _d.giveVec() + 2 * _dim, this->_mesh.numVertices() );
wr_medit_ascii2( namedef, this->_mesh, _d, this->_factor );
// wr_medit_ascii_vector("veloc."+name+".bb",_u.giveVec(),this->_mesh.numVertices(),_dim_u);
wr_medit_ascii_scalar( "dep_x." + name + ".bb", _d.giveVec(), this->mesh().numVertices() );
wr_medit_ascii_scalar( "dep_y." + name + ".bb", _d.giveVec() + _dim, this->mesh().numVertices() );
wr_medit_ascii_scalar( "dep_z." + name + ".bb", _d.giveVec() + 2 * _dim, this->mesh().numVertices() );
wr_medit_ascii2( namedef, this->mesh(), _d, this->_factor );
// wr_medit_ascii_vector("veloc."+name+".bb",_u.giveVec(),this->mesh().numVertices(),_dim_u);
system( ( "ln -s " + namedef + " dep_x." + name + ".mesh" ).data() );
system( ( "ln -s " + namedef + " dep_y." + name + ".mesh" ).data() );
system( ( "ln -s " + namedef + " dep_z." + name + ".mesh" ).data() );
// system(("ln -s "+this->_mesh_file+" veloc."+name+".mesh").data());
// system(("ln -s "+this->mesh()_file+" veloc."+name+".mesh").data());
wr_medit_ascii_scalar( "veld_x." + name + ".bb", _w.giveVec(), this->_mesh.numVertices() );
wr_medit_ascii_scalar( "veld_y." + name + ".bb", _w.giveVec() + _dim, this->_mesh.numVertices() );
wr_medit_ascii_scalar( "veld_z." + name + ".bb", _w.giveVec() + 2 * _dim, this->_mesh.numVertices() );
wr_medit_ascii_scalar( "veld_x." + name + ".bb", _w.giveVec(), this->mesh().numVertices() );
wr_medit_ascii_scalar( "veld_y." + name + ".bb", _w.giveVec() + _dim, this->mesh().numVertices() );
wr_medit_ascii_scalar( "veld_z." + name + ".bb", _w.giveVec() + 2 * _dim, this->mesh().numVertices() );
// // wr_medit_ascii_vector("veloc."+name+".bb",_u.giveVec(),this->_mesh.numVertices(),_dim_u);
// // wr_medit_ascii_vector("veloc."+name+".bb",_u.giveVec(),this->mesh().numVertices(),_dim_u);
system( ( "ln -s " + namedef + " veld_x." + name + ".mesh" ).data() );
system( ( "ln -s " + namedef + " veld_y." + name + ".mesh" ).data() );
system( ( "ln -s " + namedef + " veld_z." + name + ".mesh" ).data() );
// system(("ln -s "+this->_mesh_file+" veloc."+name+".mesh").data());
// system(("ln -s "+this->mesh()_file+" veloc."+name+".mesh").data());
}
}
......@@ -343,10 +343,10 @@ ElasticStructureHandler<Mesh>::initialize( const Function& d0, const Function& w
ID lDof;
// Loop on elements of the mesh
for ( ID iElem = 1; iElem <= this->_mesh.numVolumes(); ++iElem )
for ( ID iElem = 1; iElem <= this->mesh().numVolumes(); ++iElem )
{
_fe.updateJac( this->_mesh.volume( iElem ) );
_fe.updateJac( this->mesh().volume( iElem ) );
// Vertex based Dof
if ( nDofpV )
......@@ -455,12 +455,12 @@ ElasticStructureHandler<Mesh>::initialize( const std::string& depName,
{
std::cout << " S- restarting at time = " << startT << std::endl;
_count = (int) (startT/this->_dt - 0.5);
_count = (int) (startT/this->timestep() - 0.5);
// Loop on elements of the mesh
for ( ID iElem = 1; iElem <= this->_mesh.numVolumes(); ++iElem )
for ( ID iElem = 1; iElem <= this->mesh().numVolumes(); ++iElem )
{
_fe.updateJac( this->_mesh.volume( iElem ) );
_fe.updateJac( this->mesh().volume( iElem ) );
}
readUnknown(depName, _d);
readUnknown(velName, _w);
......
......@@ -59,9 +59,9 @@ class VenantKirchhofSolver:
{
public:
typedef ElasticStructureHandler<Mesh> super;
typedef typename super::Function Function;
typedef typename super::source_type source_type;
typedef ElasticStructureHandler<Mesh> super;
typedef typename super::Function Function;
typedef typename super::source_type source_type;
typedef typename super::bchandler_type bchandler_type;
//! Constructors
......@@ -266,13 +266,13 @@ VenantKirchhofSolver( const GetPot& data_file, const RefFE& refFE, const QuadRul
UInt nc = this->_d.nbcomp();
//inverse of dt:
Real dti2 = 2.0 / ( this->_dt * this->_dt );
Real dti2 = 2.0 / ( this->timestep() * this->timestep() );
// Elementary computation and matrix assembling
// Loop on elements
for ( UInt i = 1; i <= this->_mesh.numVolumes(); i++ )
for ( UInt i = 1; i <= this->mesh().numVolumes(); i++ )
{
this->_fe.updateFirstDerivQuadPt( this->_mesh.volumeList( i ) );
this->_fe.updateFirstDerivQuadPt( this->mesh().volumeList( i ) );
_elmatK.zero();
_elmatM.zero();
......@@ -359,14 +359,14 @@ VenantKirchhofSolver( const GetPot& data_file,
UInt nc = this->_d.nbcomp();
//inverse of dt:
Real dti2 = 2.0 / ( this->_dt * this->_dt );
Real dti2 = 2.0 / ( this->timestep() * this->timestep() );
// Elementary computation and matrix assembling
// Loop on elements
for ( UInt i = 1; i <= this->_mesh.numVolumes(); i++ )
for ( UInt i = 1; i <= this->mesh().numVolumes(); i++ )
{
this->_fe.updateFirstDerivQuadPt( this->_mesh.volumeList( i ) );
this->_fe.updateFirstDerivQuadPt( this->mesh().volumeList( i ) );
_elmatK.zero();
_elmatM.zero();
......@@ -430,10 +430,10 @@ timeAdvance( source_type const& source, const Real& time )
{
// l`oop on volumes: assembling source term
for ( UInt i = 1; i <= this->_mesh.numVolumes(); ++i )
for ( UInt i = 1; i <= this->mesh().numVolumes(); ++i )
{
this->_fe.updateFirstDerivQuadPt( this->_mesh.volumeList( i ) );
this->_fe.updateFirstDerivQuadPt( this->mesh().volumeList( i ) );
_elmatK.zero();
......@@ -466,10 +466,10 @@ timeAdvance( source_type const& source, const Real& time )
_rhsWithoutBC = ZeroVector( _rhsWithoutBC.size() );
// loop on volumes: assembling source term
for ( UInt i = 1; i <= this->_mesh.numVolumes(); ++i )
for ( UInt i = 1; i <= this->mesh().numVolumes(); ++i )
{
this->_fe.updateFirstDerivQuadPt( this->_mesh.volumeList( i ) );
this->_fe.updateFirstDerivQuadPt( this->mesh().volumeList( i ) );
_elvec.zero();
......@@ -481,11 +481,11 @@ timeAdvance( source_type const& source, const Real& time )
}
// right hand side without boundary load terms
Vector __z = this->_d + this->_dt * this->_w;
Vector __z = this->_d + this->timestep() * this->_w;
_rhsWithoutBC += _M * __z;
_rhsWithoutBC -= _K * this->_d;
_rhs_w = ( 2.0 / this->_dt ) * this->_d + this->_w;
_rhs_w = ( 2.0 / this->timestep() ) * this->_d + this->_w;
std::cout << std::endl;
std::cout << "rhsWithoutBC norm = " << norm_2(_rhsWithoutBC) << std::endl;
std::cout << "_rhs_w norm = " << norm_2(_rhs_w) << std::endl;
......@@ -517,17 +517,13 @@ iterate()
}
else
{
std::cout << "Number of inner iterations : " << maxiter << std::endl;
std::cout << "Number of inner iterations : " << maxiter << " " << this->timestep() << std::endl;
_out_iter << _time << " " << maxiter << std::endl;
}
this->_w = ( 2.0 / this->_dt ) * this->_d - _rhs_w;
this->_w = ( 2.0 / this->timestep() ) * this->_d - _rhs_w;
// evalResidual(_residual_d, _d, 0);
// _residual_d = -1*(_C*this->_d);
// std::cout << "rhsWithoutBC norm = " << norm_2(_rhsWithoutBC) << std::endl;
_residual_d = _C*this->_d - _rhsWithoutBC;
// _residual_d = -1.*_residual_d;
std::cout << " ok. " << std::flush;
}
......@@ -555,7 +551,7 @@ iterate(Vector &_sol)
_out_iter << _time << " " << maxiter << std::endl;
}
this->_w = ( 2.0 / this->_dt ) * this->_d - _rhs_w;
this->_w = ( 2.0 / this->timestep() ) * this->_d - _rhs_w;
std::cout << "sol norm = " << norm(this->_sol) << std::endl;
......@@ -593,10 +589,10 @@ evalResidual( Vector &res, const Vector& sol, int /*iter*/)
// Elementary computation and matrix assembling
// Loop on elements
for ( UInt i = 1; i <= this->_mesh.numVolumes(); i++ )
for ( UInt i = 1; i <= this->mesh().numVolumes(); i++ )
{
this->_fe.updateFirstDerivQuadPt( this->_mesh.volumeList( i ) );
this->_fe.updateFirstDerivQuadPt( this->mesh().volumeList( i ) );
_elmatK.zero();
......@@ -626,14 +622,14 @@ evalResidual( Vector &res, const Vector& sol, int /*iter*/)
std::cout << "updating the boundary conditions" << std::flush;
if ( !this->BCh_solid().bdUpdateDone() )
this->BCh_solid().bdUpdate( this->_mesh, this->_feBd, this->_dof );
this->BCh_solid().bdUpdate( this->mesh(), this->_feBd, this->_dof );
std::cout << std::endl;
bcManageMatrix( _K, this->_mesh, this->_dof, this->BCh_solid(), this->_feBd, 1.0 );
bcManageMatrix( _K, this->mesh(), this->_dof, this->BCh_solid(), this->_feBd, 1.0 );
_rhs = _rhsWithoutBC;
bcManageVector( _rhs, this->_mesh, this->_dof, this->BCh_solid(), this->_feBd, _time, 1.0 );
bcManageVector( _rhs, this->mesh(), this->_dof, this->BCh_solid(), this->_feBd, _time, 1.0 );
res = _K * sol - _rhs;
......@@ -666,9 +662,9 @@ updateJacobian( Vector& sol, int iter )
UInt nc = this->_d.nbcomp();
// loop on volumes: assembling source term
for ( UInt i = 1; i <= this->_mesh.numVolumes(); ++i )
for ( UInt i = 1; i <= this->mesh().numVolumes(); ++i )
{
this->_fe.updateFirstDerivQuadPt( this->_mesh.volumeList( i ) );
this->_fe.updateFirstDerivQuadPt( this->mesh().volumeList( i ) );
_elmatK.zero();
......@@ -723,9 +719,9 @@ solveJac( Vector &step, const Vector& res, double& /*linear_rel_tol*/)
// BC manage for the velocity
if ( !this->BCh_solid().bdUpdateDone() )
this->BCh_solid().bdUpdate( this->_mesh, this->_feBd, this->_dof );
this->BCh_solid().bdUpdate( this->mesh(), this->_feBd, this->_dof );
bcManageMatrix( _J, this->_mesh, this->_dof, this->BCh_solid(), this->_feBd, tgv );
bcManageMatrix( _J, this->mesh(), this->_dof, this->BCh_solid(), this->_feBd, tgv );
chrono.stop();
std::cout << "done in " << chrono.diff() << "s." << std::endl;
......@@ -764,9 +760,9 @@ solveJac( Vector &step, const Vector& res, double& /*linear_rel_tol*/,
// BC manage for the velocity
if ( BCh->bdUpdateDone() )
BCh->bdUpdate( this->_mesh, this->_feBd, this->_dof );
BCh->bdUpdate( this->mesh(), this->_feBd, this->_dof );
bcManageMatrix( _J, this->_mesh, this->_dof, *BCh, this->_feBd, tgv );
bcManageMatrix( _J, this->mesh(), this->_dof, *BCh, this->_feBd, tgv );
chrono.stop();
std::cout << "done in " << chrono.diff() << "s." << std::endl;
......@@ -803,9 +799,9 @@ solveJac( Vector &step, const Vector& res, double& /*linear_rel_tol*/,
// // BC manage for the velocity
// if ( !BCd.bdUpdateDone() )
// BCd.bdUpdate( this->_mesh, this->_feBd, this->_dof );
// BCd.bdUpdate( this->mesh(), this->_feBd, this->_dof );
// bcManageMatrix( _J, this->_mesh, this->_dof, BCd, this->_feBd, tgv );
// bcManageMatrix( _J, this->mesh(), this->_dof, BCd, this->_feBd, tgv );
// chrono.stop();
// std::cout << "done in " << chrono.diff() << "s." << std::endl;
......@@ -847,16 +843,16 @@ solveJacobian( Real /*time*/ )
// BC manage for the velocity
if ( !this->BCh_solid().bdUpdateDone() )
this->BCh_solid().bdUpdate( this->_mesh, this->_feBd, this->_dof );
this->BCh_solid().bdUpdate( this->mesh(), this->_feBd, this->_dof );
bcManageVector(_f,
this->_mesh,
this->mesh(),
this->_dof,
this->BCh_solid(),
this->_feBd,
1., 1.);
bcManageMatrix( _J, this->_mesh, this->_dof, this->BCh_solid(), this->_feBd, tgv );
bcManageMatrix( _J, this->mesh(), this->_dof, this->BCh_solid(), this->_feBd, tgv );
chrono.stop();
std::cout << "done in " << chrono.diff() << "s." << std::endl;
......@@ -869,7 +865,7 @@ solveJacobian( Real /*time*/ )
chrono.stop();
std::cout << "done in " << chrono.diff() << " s." << std::endl;
this->_w = ( 2.0 / this->_dt ) * M_ddisp - _rhs_w;
this->_w = ( 2.0 / this->timestep() ) * M_ddisp - _rhs_w;
// std::cout << " S- Computing residual ... " << std::flush;
_residual_d = _C*M_ddisp;
......@@ -896,16 +892,16 @@ solveJacobian( const Real /*time*/ , bchandler_type& BCd)
// BC manage for the velocity
if ( !(*BCd).bdUpdateDone() )
(*BCd).bdUpdate( this->_mesh, this->_feBd, this->_dof );
(*BCd).bdUpdate( this->mesh(), this->_feBd, this->_dof );
bcManageVector(_f,
this->_mesh,
this->mesh(),
this->_dof,
*BCd,
this->_feBd,
1., 1.);
bcManageMatrix( _J, this->_mesh, this->_dof, *BCd, this->_feBd, tgv );
bcManageMatrix( _J, this->mesh(), this->_dof, *BCd, this->_feBd, tgv );
chrono.stop();
std::cout << "done in " << chrono.diff() << "s." << std::endl;
......@@ -918,7 +914,7 @@ solveJacobian( const Real /*time*/ , bchandler_type& BCd)
chrono.stop();
std::cout << "done in " << chrono.diff() << " s." << std::endl;
this->_w = ( 2.0 / this->_dt ) * M_ddisp - _rhs_w;
this->_w = ( 2.0 / this->timestep() ) * M_ddisp - _rhs_w;
// std::cout << " S- Computing residual ... " << std::flush;
_residual_d = _C*M_ddisp;
......
......@@ -170,16 +170,16 @@ ConvDiffReactHandler( const GetPot& data_file, const RefFE& refFE_c,
const QuadRule& Qr_c, const QuadRule& bdQr_c, BCHandler& BCh_c ) :
DataConvDiffReact<Mesh>( data_file ),
_refFE_c( refFE_c ),
_dof_c( this->_mesh, _refFE_c ),
_dof_c( this->mesh(), _refFE_c ),
_dim_c( _dof_c.numTotalDof() ),
_Qr_c( Qr_c ),
_bdQr_c( bdQr_c ),
_fe_c( _refFE_c, getGeoMap( this->_mesh ), _Qr_c ),
_feBd_c( _refFE_c.boundaryFE(), getGeoMap( this->_mesh ).boundaryMap(), _bdQr_c ),
_fe_c( _refFE_c, getGeoMap( this->mesh() ), _Qr_c ),
_feBd_c( _refFE_c.boundaryFE(), getGeoMap( this->mesh() ).boundaryMap(), _bdQr_c ),
_c( _dim_c ),
_u_c( _dim_c ),
_BCh_c( BCh_c ),
_bdf( this->_order_bdf )
_bdf( this->order_bdf() )
{}
......@@ -223,7 +223,7 @@ void
ConvDiffReactHandler<Mesh>::initialize( const Function& c0, Real t0, Real dt )
{
_bdf.initialize_unk( c0, this->_mesh, _refFE_c, _fe_c, _dof_c, t0, dt, 1 );
_bdf.initialize_unk( c0, this->mesh(), _refFE_c, _fe_c, _dof_c, t0, dt, 1 );
_c = *( _bdf.unk().begin() );
_bdf.showMe();
......
......@@ -168,7 +168,7 @@ ConvDiffReactSolverPC( const GetPot& data_file, const RefFE& refFE_c, const Quad
_M_c.zeros();
//inverse of dt:
Real dti = 1. / this->_dt;
Real dti = 1. / this->timestep();
// *******************************************************
// Coefficient of the mass term at time t^{n+1}
......@@ -180,10 +180,10 @@ ConvDiffReactSolverPC( const GetPot& data_file, const RefFE& refFE_c, const Quad
// Elementary computation and matrix assembling
for ( UInt i = 1; i <= this->_mesh.numVolumes(); i++ )
for ( UInt i = 1; i <= this->mesh().numVolumes(); i++ )
{ // Loop on elements
this->_fe_c.updateFirstDerivQuadPt( this->_mesh.volumeList( i ) );
this->_fe_c.updateFirstDerivQuadPt( this->mesh().volumeList( i ) );
_elmatC.zero();
_elmatM_c.zero();
......@@ -230,10 +230,10 @@ timeAdvance( source_type const& source, const Real& time )
_f_c = ZeroVector( _f_c.size() );
// loop on volumes: assembling source term
for ( UInt i = 1; i <= this->_mesh.numVolumes(); ++i )
for ( UInt i = 1; i <= this->mesh().numVolumes(); ++i )
{
_elvec.zero();
this->_fe_c.update( this->_mesh.volumeList( i ) );
this->_fe_c.update( this->mesh().volumeList( i ) );
compute_vec( source, _elvec, this->_fe_c, time, 0 ); // compute local vector
assemb_vec( _f_c, _elvec, this->_fe_c, this->_dof_c, 0 ); // assemble local vector into global one
......@@ -269,10 +269,10 @@ iterate( const Real& time )
chrono.start();
// loop on volumes
for ( UInt i = 1; i <= this->_mesh.numVolumes(); ++i )
for ( UInt i = 1; i <= this->mesh().numVolumes(); ++i )
{
this->_fe_c.updateFirstDeriv( this->_mesh.volumeList( i ) ); // as updateFirstDer
this->_fe_c.updateFirstDeriv( this->mesh().volumeList( i ) ); // as updateFirstDer
_elmatC.zero();
......@@ -352,8 +352,8 @@ iterate( const Real& time )
chrono.start();
// BC manage for the concentration
if ( !this->_BCh_c.bdUpdateDone() )
this->_BCh_c.bdUpdate( this->_mesh, this->_feBd_c, this->_dof_c );
bcManage( _CDR, _f_c, this->_mesh, this->_dof_c, this->_BCh_c, this->_feBd_c, tgv, time );
this->_BCh_c.bdUpdate( this->mesh(), this->_feBd_c, this->_dof_c );
bcManage( _CDR, _f_c, this->mesh(), this->_dof_c, this->_BCh_c, this->_feBd_c, tgv, time );
chrono.stop();
std::cout << "done in " << chrono.diff() << "s." << std::endl;
......@@ -409,15 +409,15 @@ getvel( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u, con
{
for ( UInt i = 0; i < this->_mesh.numVertices(); i++ )
for ( UInt i = 0; i < this->mesh().numVertices(); i++ )
{
if ( this->_u_to_c[ i ].ele == 0 )
{
// Dirichlet boundary for the velocity -> get velocity for boundary function
Real xp = this->_mesh.point( i + 1 ).x();
Real yp = this->_mesh.point( i + 1 ).y();
Real zp = this->_mesh.point( i + 1 ).z();
Real xp = this->mesh().point( i + 1 ).x();
Real yp = this->mesh().point( i + 1 ).y();
Real zp = this->mesh().point( i + 1 ).z();
for ( ID jj = 0; jj < 3; ++jj )
this->_u_c( i + jj * this->_u_c.size() / 3 ) = BCh_u[ ( int ) this->_u_to_c[ i ].b[ 0 ] ] ( time, xp, yp, zp, BCh_u[ ( int ) this->_u_to_c[ i ].b[ 0 ] ].component( jj + 1 ) );
......@@ -470,17 +470,17 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
SimpleVect<GeoElement3D<LinearTetra> >::iterator iv = umesh.volumeList.begin();
for ( UInt i = 0; i < this->_mesh.numVertices(); i++ )
for ( UInt i = 0; i < this->mesh().numVertices(); i++ )
{
tetrapoints( 1 ) = this->_mesh.point( i + 1 );
tetrapoints( 1 ) = this->mesh().point( i + 1 );
if ( this->_mesh.point( i + 1 ).boundary() )
if ( this->mesh().point( i + 1 ).boundary() )
{
UInt l;
for ( UInt k = 0; k < BCh_u.size(); k++ )
{
if ( BCh_u[ k ].flag() == this->_mesh.point( i + 1 ).marker() )
if ( BCh_u[ k ].flag() == this->mesh().point( i + 1 ).marker() )
{
l = k;
}
......@@ -525,7 +525,7 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
}
if (volume < minvolume)
{
if (adjacent != 0)
if (adjacent != 0)
{
minvolume = volume;
minadjacent=adjacent;
......@@ -551,7 +551,6 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
localcoord.b[ 3 ] = b3;
localcoord.ele = vid;
this->_u_to_c.push_back( localcoord );
}
else
{
......@@ -563,14 +562,12 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
tetrapoints_found( 4 ) = umesh.point( ( iv->point( 4 ) ).id() );
test( tetrapoints_found, tetrapoints, b1, b2, b3 );
localcoord.b[ 0 ] = 1.0 - b1 - b2 - b3;
localcoord.b[ 1 ] = b1;
localcoord.b[ 2 ] = b2;
localcoord.b[ 3 ] = b3;
localcoord.ele = vid;
this->_u_to_c.push_back( localcoord );
}
}
}
......@@ -608,7 +605,7 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
}
if (volume < minvolume)
{
if (adjacent != 0)
if (adjacent != 0)
{
minvolume = volume;
minadjacent=adjacent;
......@@ -636,7 +633,6 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
localcoord.b[ 3 ] = b3;
localcoord.ele = vid;
this->_u_to_c.push_back( localcoord );
}
else
{
......@@ -647,14 +643,12 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
tetrapoints_found( 4 ) = umesh.point( ( iv->point( 4 ) ).id() );
test( tetrapoints_found, tetrapoints, b1, b2, b3 );
localcoord.b[ 0 ] = 1.0 - b1 - b2 - b3;
localcoord.b[ 1 ] = b1;
localcoord.b[ 2 ] = b2;
localcoord.b[ 3 ] = b3;
localcoord.ele = vid;
this->_u_to_c.push_back( localcoord );
}
}
}
......
......@@ -148,7 +148,7 @@ DarcyHandler<Mesh>::DarcyHandler( const GetPot& data_file, const RefHdivFE& refF
DataDarcy<Mesh>(data_file),
DataAztec(data_file, "darcy/aztec"),
nbCoor(nDimensions),
geoMap( getGeoMap( this->_mesh ) ),
geoMap( getGeoMap( this->mesh() ) ),
qr( qr_u ),
geoMapBd( geoMap.boundaryMap() ),
qrBd( bdqr_u ),
......@@ -164,20 +164,20 @@ DarcyHandler<Mesh>::DarcyHandler( const GetPot& data_file, const RefHdivFE& refF
vdof( refVFE ),
pdof( refPFE ),
tpdof( refTPFE ),
numFacesPerVolume( this->_mesh.volumeList(1).numLocalFaces )
numFacesPerVolume( this->mesh().volumeList(1).numLocalFaces )
/* we assume that all element have the same number
of faces, so we just look at the first one */
{
//! This function should be already called in DataMesh construction
// this->_mesh.updateElementFaces(true);
// this->mesh().updateElementFaces(true);
/*the "true" flag is to build the faceList
of all faces (and not only the boundary) */
if(this->verbose>2) this->_mesh.showMe();
if(this->verbose>2) this->mesh().showMe();
// build dof
vdof.update(this->_mesh);
pdof.update(this->_mesh);
tpdof.update(this->_mesh);
vdof.update(this->mesh());
pdof.update(this->mesh());
tpdof.update(this->mesh());
dimTPdof = tpdof.numTotalDof();
dimPdof = pdof.numTotalDof();
......@@ -192,7 +192,7 @@ DarcyHandler<Mesh>::DarcyHandler( const GetPot& data_file, const RefHdivFE& refF
/*
// check the mesh after b.c.
BE CAREFUL: calling
this->_mesh.check(true,true);
this->mesh().check(true,true);
after updateElementFaces(true) erase faceElement !!!!
*/
//
......
......@@ -112,7 +112,7 @@ public:
_M_bc = __bch;
// update the dof with the b.c.
_M_bc.bdUpdate(this->_mesh, this->feBd, this->tpdof);