Commit 0158b2db authored by mm0's avatar mm0
Browse files

some little bug fixes and added in testsuite/lifefilters/ tests about netgen

parent 9083e62e
......@@ -5,6 +5,10 @@ Bugs fixed in 0.5.0:
[configure]
* Fixed umfpack check when umfpack is not installed (CP)
[all]
* fixed little bugs in convDiffReactSolverPC,
basisElSh,currentFE,defQuadRuleFE (MM)
New in 0.5.0:
=============
[configure]
......@@ -44,6 +48,7 @@ New in 0.5.0:
* Removed test_darcytetra (VM)
* oned solver can solve 2 coupled tubes (VM)
* Centralized all FSI tests in test_fsi (GF)
* Added test for netgen in lifefilters (MM)
Bugs fixed in 0.4.0:
====================
......
......@@ -134,9 +134,9 @@ typedef size_t ID;
//typedef unsigned int ID;
//! I define UInt=unsigned int. This allow to use this stuff also outside lifeV
// more correct version
typedef size_t UInt;
//typedef size_t UInt;
// original version
//typedef unsigned int UInt;
typedef unsigned int UInt;
#endif
//! The Edge basis class
......
......@@ -20,7 +20,7 @@ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
/*! Contains the basic element shapes, to be used by Geometric and Finite
Elements
$Header: /cvsroot/lifev/lifev/life/lifecore/Attic/basisElSh.hpp,v 1.6 2004-11-08 10:16:19 prudhomm Exp $
$Header: /cvsroot/lifev/lifev/life/lifecore/Attic/basisElSh.hpp,v 1.7 2004-11-14 18:11:45 mm0 Exp $
\version 0.0 Experimental 19/8/99. Luca Formaggia
......@@ -104,7 +104,7 @@ public:
static const ReferenceGeometry Geometry = EDGE;
static const UInt nDim = 1;
static const UInt numFaces = 0;
static const UInt numEdges = 0;
static const UInt numEdges = 1;
static const UInt numVertices = 2;
};
......@@ -116,7 +116,7 @@ public:
static const ReferenceGeometry Geometry = FACE;
static const UInt nDim = 2;
static const UInt numVertices = 3;
static const UInt numFaces = 0;
static const UInt numFaces = 1;
static const UInt numEdges = numVertices;
};
......@@ -127,7 +127,7 @@ public:
static const ReferenceShapes Shape = QUAD;
static const ReferenceGeometry Geometry = FACE;
static const UInt nDim = 2;
static const UInt numFaces = 0;
static const UInt numFaces = 1;
static const UInt numVertices = 4;
static const UInt numEdges = numVertices;
};
......
......@@ -511,6 +511,7 @@ void CurrentFE::_comp_phiDerDer2()
x1 = 0.;
for ( int jcoor = 0;jcoor < nbCoor;jcoor++ )
{
x1 += tInvJac(icoor,jcoor,ig)*dPhiRef(j,jcoor,ig);
x2 = 0.;
for ( int k1 = 0;k1 < nbCoor;k1++ )
{
......
......@@ -529,11 +529,11 @@ Real fct1_P2_1D( cRRef x, cRRef, cRRef )
{
return 2. * ( x - 1. ) * ( x - 0.5 );
}
Real fct2_P2_1D( cRRef x, cRRef, cRRef )
Real fct3_P2_1D( cRRef x, cRRef, cRRef )
{
return 4. * x * ( 1. - x );
}
Real fct3_P2_1D( cRRef x, cRRef, cRRef )
Real fct2_P2_1D( cRRef x, cRRef, cRRef )
{
return 2. * x * ( x - 0.5 );
}
......@@ -542,11 +542,11 @@ Real derfct1_1_P2_1D( cRRef x, cRRef , cRRef )
{
return 4. * x - 3.;
}
Real derfct2_1_P2_1D( cRRef x, cRRef , cRRef )
Real derfct3_1_P2_1D( cRRef x, cRRef , cRRef )
{
return -8. * x + 4.;
}
Real derfct3_1_P2_1D( cRRef x, cRRef , cRRef )
Real derfct2_1_P2_1D( cRRef x, cRRef , cRRef )
{
return 4. * x - 1.;
}
......@@ -555,11 +555,11 @@ Real der2fct1_11_P2_1D( cRRef x, cRRef , cRRef )
{
return 4;
}
Real der2fct2_11_P2_1D( cRRef x, cRRef , cRRef )
Real der2fct3_11_P2_1D( cRRef x, cRRef , cRRef )
{
return -8;
}
Real der2fct3_11_P2_1D( cRRef x, cRRef , cRRef )
Real der2fct2_11_P2_1D( cRRef x, cRRef , cRRef )
{
return 4;
}
......
......@@ -364,7 +364,7 @@ bool checkMesh3D( RegionMesh3D & mesh, Switch & sw,
for ( typename RegionMesh3D::Volumes::iterator iv = mesh.volumeList.begin();
iv != mesh.volumeList.end();++iv )
{
if ( iv->markerUnset() )
if ( iv->isMarkerUnset() )
iv->setMarker( mesh.marker() );
}
}
......
This diff is collapsed.
......@@ -75,13 +75,14 @@ Geo0D::showMe( bool verbose, std::ostream & out ) const
out << " Geo0D object " << std::endl;
if ( verbose )
{
unsigned i;
out << " Coordinates:" << std::endl;
Real const * c = coor();
for ( unsigned i = 0; i < nDimensions; i++ )
for ( i = 0; i < nDimensions-1; i++ )
{
out << c[ i ] << ", ";
}
out << std::endl << std::endl;
out << c[i] << std::endl << std::endl;
}
out << "ID= " << id() << " ";
out << "----- END OF Geo0D data ---" << std::endl << std::endl;
......
......@@ -27,17 +27,18 @@ namespace LifeV
//MarkerTraits_Base
const MarkerTraits_Base::EntityFlag MarkerTraits_Base::NULLFLAG = LONG_MIN;
//MM: if you modity these changes here recheck function readNetgenMesh
// becouse it uses this changes
MarkerTraits_Base::EntityFlag MarkerTraits_Base::strongerFlag( EntityFlag const & a, EntityFlag const & b )
{
if ( a == NULLFLAG | b == NULLFLAG )
return NULLFLAG;
return a > b ? a : b ;
}
MarkerTraits_Base::EntityFlag MarkerTraits_Base::weakerFlag( EntityFlag const & a, EntityFlag const & b )
{
if ( a == NULLFLAG | b == NULLFLAG )
return NULLFLAG;
if(a==NULLFLAG)return b;
if(b==NULLFLAG)return a;
return a < b ? a : b ;
}
}
......@@ -166,10 +166,10 @@ public:
inline bool isMarkerUnset() const;
//! Put marker to nullflag
inline bool unsetMarker() const;
inline void unsetMarker() const;
//! Put marker to nullflag
inline bool markerUnset() const;
inline void markerUnset() const;
//! Helper function that prints a marker Flag
std::ostream & printFlag( EntityFlag const f, std::ostream & out ) const;
......@@ -284,7 +284,7 @@ typename MarkerTraits::EntityFlag Marker_Base<MarkerTraits>::setWeakerMarker( En
template <typename MarkerTraits>
typename MarkerTraits::EntityFlag Marker_Base<MarkerTraits>::setStrongerMarker( EntityFlag const & p )
{
if ( markerUnset() )
if ( isMarkerUnset() )
return flag = p;
return setMarker( MarkerTraits::strongerFlag( this->marker(), p ) );
}
......@@ -292,7 +292,7 @@ typename MarkerTraits::EntityFlag Marker_Base<MarkerTraits>::setStrongerMarker(
template <typename MarkerTraits>
typename MarkerTraits::EntityFlag Marker_Base<MarkerTraits>::setWeakerMarker( EntityFlag const & p )
{
if ( markerUnset() )
if ( isMarkerUnset() )
return flag = p;
return setMarker( MarkerTraits::weakerFlag( this->marker(), p ) );
}
......@@ -310,16 +310,16 @@ bool Marker_Base<MarkerTraits>::isMarkerUnset() const
}
template <typename MarkerTraits>
bool Marker_Base<MarkerTraits>::unsetMarker() const
void Marker_Base<MarkerTraits>::unsetMarker() const
{
return isMarkerUnset();
flag=nullFlag();
}
template <typename MarkerTraits>
bool Marker_Base<MarkerTraits>::markerUnset() const
void Marker_Base<MarkerTraits>::markerUnset() const
{
return isMarkerUnset();
flag=nullFlag();
}
template <typename MarkerTraits>
......
......@@ -1025,7 +1025,7 @@ bool fixBoundaryFaces( RegionMesh3D & mesh,
// Correct extra info
fit->ad_first() = vol;
fit->pos_first() = j;
if ( fit->markerUnset() )
if ( fit->isMarkerUnset() )
{
inheritWeakerMarker( *fit );
if ( verbose )
......@@ -1532,7 +1532,7 @@ p1top2( RegionMesh & mesh, std::ostream & out = std::cout )
//#define JFG
//#ifndef JFG
// original version: DOES NOT work when the edges are not give in the mesh file
if ( pe->markerUnset() )
if ( pe->isMarkerUnset() )
inheritWeakerMarker( *pe );
pp->setMarker( pe->marker() );
//#else
......
......@@ -524,10 +524,10 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
z[ 1 ] = umesh.point( i1 ).z();
z[ 2 ] = umesh.point( i2 ).z();
z[ 3 ] = umesh.point( i3 ).z();
volume[ jj ] = calcvol( x, y, z );
if ( volume[ jj ] < minvolume )
volume[ jj-1 ] = calcvol( x, y, z );
if ( volume[ jj-1 ] < minvolume )
{
minvolume = volume[ jj ];
minvolume = volume[ jj-1 ];
jk = jj;
}
}
......@@ -611,10 +611,10 @@ getcoord( RegionMesh3D & umesh, PhysVectUnknown<Vector> & u, BCHandler& BCh_u )
z[ 1 ] = umesh.point( i1 ).z();
z[ 2 ] = umesh.point( i2 ).z();
z[ 3 ] = umesh.point( i3 ).z();
volume[ jj ] = calcvol( x, y, z );
if ( volume[ jj ] < minvolume )
volume[ jj-1 ] = calcvol( x, y, z );
if ( volume[ jj-1 ] < minvolume )
{
minvolume = volume[ jj ];
minvolume = volume[ jj-1 ];
jk = jj;
}
}
......
......@@ -33,7 +33,8 @@ SUBDIRS = data \
test_ns_bdf test_ns_cyl test_ns_cyl_flux test_ns_bcvec test_ns_sstress \
test_fsi \
test_masstransport test_mass_coupling test_ns_ipstab \
test_coupling_3d0d
test_coupling_3d0d \
lifefilters
TESTSUITES=$(srcdir)/lifecore/testsuite.at \
$(srcdir)/lifearray/testsuite.at \
......@@ -61,7 +62,8 @@ $(srcdir)/test_ns_sstress/testsuite.at \
$(srcdir)/test_onedmodel/testsuite.at \
$(srcdir)/test_p2/testsuite.at \
$(srcdir)/test_postproc/testsuite.at \
$(srcdir)/test_q1/testsuite.at
$(srcdir)/test_q1/testsuite.at \
$(srcdir)/lifefilters/testsuite.at
$(srcdir)/package.m4: $(top_srcdir)/configure.in.in
......
This diff is collapsed.
# -*- makefile -*-
#
#
#
# Copyright (C) 2004 Christophe Prud'homme (christophe.prudhomme@epfl.ch)
#
# Distributed under the GPL(GNU Public License):
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
SUFFIXES = .cpp .hpp .idl .c .h .f .F .o .moc
noinst_HEADERS = main.hpp ud_functions.hpp
check_PROGRAMS = test_netgen
test_netgen_SOURCES = main.cpp
EXTRA_DIST = data testsuite.at
MOSTLYCLEANFILES=scal.vtk
include $(top_srcdir)/testsuite/Makefile.testsuite
# -*- getpot -*- (GetPot mode activation for emacs)
#-------------------------------------------------
# Data file for Mesh test
#-------------------------------------------------
mesh_type = 'NETGEN'
mesh_dir = '../data/mesh/netgen/' # the directory where the mesh file is
mesh_file = 'cubeBc312_x4.vol' # mesh file
#mesh_file = 'cubeBc312_x16.vol' # mesh file
mesh_faces = all # update all faces elements
/* -*- mode: c++ -*-
This program is part of the LifeV library
Copyright (C) 2001,2002,2003,2004 EPFL, INRIA, Politechnico di Milano
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/* ========================================================
Simple Laplacian test with Dirichlet Boundary condition
Solve the problem
- \Delta u = f
u = g on the boundary
on a cube
*/
/* Mauro Meneghin: note, I've taken the test_p2 example and
I've modified it a bit to read netgen meshes too,
this files recognizes old meshes too anyway, so you could
use this file as a template instead of the old ones */
#include <GetPot.hpp>
#include "main.hpp"
#include "ud_functions.hpp"
//#include "bc_manage.hpp" bc_manage.hpp changed to bcManage.hpp,
//BCFunction_Base to BCFunctionBase, BC_Handler to BCHandler
#include "bcManage.hpp"
#include "elemMat.hpp"
#include "elemOper.hpp"
#include "openDX_wrtrs.hpp"
#include "vtk_wrtrs.hpp"
//#include "sobolevNorms.hpp"
#undef OPER_TEMPLATE
//#define P1
#define P2
#undef INRIA
int main() {
using namespace LifeV;
using namespace std;
Chrono chrono;
// ===================================================
// Boundary conditions definition
// ===================================================
BCFunctionBase gv1(g1); // Functor storing the user definded function g
BCFunctionBase gv3(g3); // Functor storing the user definded function g
BCFunctionBase gv2(g2); // Functor storing the user definded function g
BCHandler BCh(3); // We impose three boundary conditions
BCh.addBC("Adiabatic", 1, Natural, Scalar,gv1);
BCh.addBC("Left", 3, Essential, Scalar, gv3);
BCh.addBC("Right", 2, Essential, Scalar, gv2);
// Ouput
BCh.showMe();
// ===================================================
// Finite element stuff
// ===================================================
const GeoMap& geoMap = geoLinearTetra;
const QuadRule& qr = quadRuleTetra5pt;
const GeoMap& geoMapBd = geoLinearTria;
const QuadRule& qrBd = quadRuleTria3pt;
#ifndef P2
// P1 elements
const RefFE& refFE = feTetraP1;
const RefFE& refBdFE = feTriaP1;
#else
//P2 elements
const RefFE& refFE = feTetraP2;
const RefFE& refBdFE = feTriaP2;
#endif
// ===================================================
// Mesh stuff
// ===================================================
#ifdef P2
// RegionMesh3D<LinearTetra> aMesh;
RegionMesh3D<QuadraticTetra> aMesh;
#else
RegionMesh3D<LinearTetra> aMesh;
#endif
long int m=1;
GetPot datafile( "data" );
std::string mesh_type = datafile( "mesh_type", "NETGEN" );
if ( mesh_type == "INRIA" )
{
string mesh_dir = datafile( "mesh_dir", "." );
string fname=mesh_dir+datafile( "mesh_file", "." );
readINRIAMeshFile(aMesh,fname,m);
}
else if ( mesh_type == "MESH++" )
{
string mesh_dir = datafile( "mesh_dir", "." );
string fname=mesh_dir+datafile( "mesh_file", "." );
// string fname=mesh_dir+"cube_48.m++";
readMppFile(aMesh,fname,m);
}
else if (mesh_type=="NETGEN"){
string mesh_dir = datafile( "mesh_dir", "." );
string fname=mesh_dir+datafile( "mesh_file", "." );
std::cout<<"opening mesh file "<<fname<<std::endl;
readNetgenMesh(aMesh,fname,m);
}
else
{
std::cerr << "wrong mesh type. It can be either MESH++ or INRIA or NETGEN" << std::endl;
return EXIT_FAILURE;
}
// Avaliable meshes
//string fname=mesh_dir+"cube_48.m++";
aMesh.showMe();
cout<< "Now building local Edges/faces Stuff"<<endl<<endl;
aMesh.updateElementEdges();
aMesh.updateElementFaces();
aMesh.showMe();
// ===================================================
// Current FE classes for the problem under study with
// mapping and quadrature rules
// ===================================================
CurrentFE fe(refFE,geoMap,qr);
CurrentBdFE feBd(refBdFE,geoMapBd,qrBd);
// ===============================================
// Update of the Dof for the particular FE problem
// and for the boundary conditions
// ===============================================
Dof dof(refFE);
dof.update(aMesh);
BCh.bdUpdate( aMesh, feBd, dof );
UInt dim = dof.numTotalDof();
dof.showMe();
// initialization of vector of unknowns and rhs
ScalUnknown<Vector> U(dim), F(dim);
U=ZeroVector( dim );
F=ZeroVector( dim );
// ==========================================
// Pattern construction and matrix assembling
// ==========================================
cout << "dim = " << dim << endl << endl;
// pattern for stiff operator
MSRPatt pattA(dof);
cout << "Values" << endl;
// A: stiff operator
MSRMatr<double> A(pattA);
cout << "*** Matrix computation : "<<endl;
chrono.start();
//
SourceFct sourceFct;
#ifdef OPER_TEMPLATE
// assembling of A: stiff operator
Stiff Ostiff(&fe);
EOStiff stiff(Ostiff);
assemble(mu*stiff+sigma*mass,aMesh,fe,dof,sourceFct,A,F);
#else
ElemMat elmat(fe.nbNode,1,1);
ElemVec elvec(fe.nbNode,1);
for(UInt i = 1; i<=aMesh.numVolumes(); i++){
fe.updateFirstDerivQuadPt(aMesh.volumeList(i));
elmat.zero();
elvec.zero();
stiff(1.,elmat,fe);
source(sourceFct,elvec,fe,0);
assemb_mat(A,elmat,fe,dof,0,0);
assemb_vec(F,elvec,fe,dof,0);
}
#endif
chrono.stop();
cout << "A has been constructed" << endl;
cout << chrono.diff() << "s." << endl;
// ====================================
// Treatment of the Boundary conditions
// ====================================
// BC manage for the velocity
cout << "*** BC Management: "<<endl;
Real tgv=1.;
chrono.start();
bcManage(A,F,aMesh,dof,BCh,feBd,tgv,0.0);
chrono.stop();
cout << chrono.diff() << "s." << endl;
// ==============================
// Reolution of the linear system
// ==============================
int proc_config[AZ_PROC_SIZE];// Processor information:
// proc_config[AZ_node] = node name
// proc_config[AZ_N_procs] = # of nodes
int options[AZ_OPTIONS_SIZE]; // Array used to select solver options.
double params[AZ_PARAMS_SIZE]; // User selected solver paramters.
int *data_org; // Array to specify data layout
double status[AZ_STATUS_SIZE]; // Information returned from AZ_solve()
// indicating success or failure.
// altre dichiarazioni per AZTEC
int *update, // vector elements updated on this node.
*external; // vector elements needed by this node.
int *update_index; // ordering of update[] and external[]
int *extern_index; // locally on this processor.
// int *bindx; // Sparse matrix to be solved is stored
// double *val; // in these MSR arrays.
int N_update; // # of unknowns updated on this node
//
cout << "*** Linear System Solving (AZTEC)" << endl;
AZ_set_proc_config(proc_config, AZ_NOT_MPI );
// cout << AZ_PROC_SIZE << " " << AZ_node << " " << AZ_N_procs << endl;
// for (UInt ii=0; ii<AZ_PROC_SIZE; ++ii)
// cout << proc_config[ii] << endl;
AZ_read_update(&N_update, &update, proc_config, U.size(), 1, AZ_linear);
AZ_defaults(options,params);
AZ_transform(proc_config, &external,
(int *)pattA.giveRaw_bindx(), A.giveRaw_value(),
update, &update_index,
&extern_index, &data_org, N_update, NULL, NULL, NULL, NULL,
AZ_MSR_MATRIX);
chrono.start();
init_options(options,params);
AZ_solve(U.giveVec(), F.giveVec(), options, params, NULL,
(int *)pattA.giveRaw_bindx(), NULL, NULL, NULL,
A.giveRaw_value(), data_org,
status, proc_config);
//
chrono.stop();