Skip to content
Snippets Groups Projects
Commit 6dc7b68a authored by Robert Klöfkorn's avatar Robert Klöfkorn
Browse files

moved to albertagrid/ 

[[Imported from SVN: r1252]]
parent 4289ea0e
Branches
Tags
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 8887 deletions
grid/albertgrid/Makefile.am deleted 100644 → 0
+ 0
− 5
View file @ 4289ea0e
# $Id$
albertdir = $(includedir)/dune/grid/albertgrid/
albert_HEADERS = agelementindex.cc albertgrid.cc part.cc partialGrid.cc \
agcommunicator.hh agindex.hh agmemory.hh albertextra.hh
grid/albertgrid/agcommunicator.hh deleted 100644 → 0
+ 0
− 578
View file @ 4289ea0e
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ALBERTGRID_COMMUNICATOR_HH
#define DUNE_ALBERTGRID_COMMUNICATOR_HH
// use this define to control if Albert should use the found MPI
// #define ALBERT_USES_MPI 1
#ifdef HAVE_CONFIG_H
#include <config.h>
#if defined(HAVE_MPI) && defined(ALBERT_USES_MPI)
#include <mpi.h>
#endif
#endif
#include <dune/common/dlist.hh>
namespace Dune {
static const int COMMUNICATOR_COMM_TAG = 457;
/*!
ProcListElement describes the link between two processors.
It contains all needed information for cummunication between
these two procs.
*/
#if defined(HAVE_MPI) && defined(ALBERT_USES_MPI)
template <class BufferType>
class ProcListElement
{
public:
//! Constructor
ProcListElement(const int owner, const int numFuncs )
: owner_ (owner) , numFuncs_ (numFuncs)
, recvCount(0), sendCount(0)
, sendBuffer_ (0) , recvBuffer_ (0)
, recive_(100) , send_(100)
{ }
//! Destructor
~ProcListElement()
{
if(sendBuffer_) std::free(sendBuffer_);sendBuffer_ = 0;
if(recvBuffer_) std::free(recvBuffer_);recvBuffer_ = 0;
}
//! return my processor number
int processor () const { return owner_; };
//! set
void setSize( int recvSize , int sendSize )
{
recvCount = recvSize;
sendCount = sendSize;
}
//! insert new recive element
void insertRecive( int num )
{
if(recvCount >= (int) recive_.size())
{
recive_.resize( 2*recvCount );
}
recive_[recvCount] = num;
recvCount++;
}
//! insert new send element
void insertSend( int num )
{
if(sendCount >= (int)send_.size())
{
send_.resize( 2*sendCount );
}
send_[sendCount] = num;
sendCount++;
}
//! write elements to buffer
template <class DofArrayType>
void loadSendBuffer (DofArrayType & vec, int pos=0)
{
int beg = pos * sendSize();
for(int i=beg; i<beg+sendSize(); i++)
{
sendBuffer_[i] = vec[send_[i-beg]];
}
}
//! read elements from buffer
template <class DofArrayType>
void unloadRecvBuffer (DofArrayType & vec, int pos=0)
{
int beg = pos * recvSize();
for(int i=beg; i<beg+recvSize(); i++)
{
vec[recive_[i-beg]] = recvBuffer_[i];
}
}
//! allocate memory for buffers
void makeBuffer()
{
recvBuffer_ = (BufferType *) std::malloc( realRecvSize()*sizeof(BufferType) );
for(int i=0 ; i<realRecvSize(); i++)
{
recvBuffer_[i] = -1.0;
}
sendBuffer_ = (BufferType *) std::malloc( realSendSize()*sizeof(BufferType) );
}
//! size of list of send elements
int sendSize () const { return sendCount; }
//! size of send buffer
int realSendSize () const { return numFuncs_ * sendSize(); }
//! pointer to send buffer
BufferType * getSendBuffer() const { return sendBuffer_; }
//! size of list of recive elements
int recvSize () const { return recvCount; }
//! size of recive buffer
int realRecvSize () const { return numFuncs_ * recvSize(); }
//! pointer to recive buffer
BufferType * getRecvBuffer() const { return recvBuffer_; }
//! print link list
void print(std::ostream &s) const
{
s << "ProcList for Processor " << owner_ << "\n";
s << "SendList " << sendSize() << "\n";
for(int i=0; i<sendSize(); i++)
{
s << send_[i] << " ";
}
s << "\n";
s << "RecvList " << recvSize() << "\n";
for(int i=0; i<recvSize(); i++)
{
s << recive_[i] << " ";
}
s << "\n";
};
//! the send and recive MPI_Request
MPI_Request * sendreq() { return &sendreq_; }
MPI_Request * recvreq() { return &recvreq_; }
private:
const int owner_;
const int numFuncs_;
int recvCount;
int sendCount;
BufferType * sendBuffer_;
BufferType * recvBuffer_;
MPI_Request sendreq_;
MPI_Request recvreq_;
// buffer size of processor
std::vector<int> recive_;
std::vector<int> send_;
};
#endif
/*!
AlbertGridCommunicator organizes the communication of AlbertGrid on
diffrent processors.
*/
template <class T>
class GathScatt
{
T fake;
public:
template <class VecType>
void gather (VecType & t, int index )
{
fake = t[index];
}
template <class VecType>
void scatter(VecType & t , int index)
{
t[index] = fake;
}
};
#if defined(HAVE_MPI) && defined(ALBERT_USES_MPI)
template <class GridType, class IndexSetType>
class AlbertGridCommunicator
{
typedef DoubleLinkedList< ProcListElement<double> > ProcessorListType;
public:
//! Constructor
AlbertGridCommunicator(GridType &grid, IndexSetType &indexSet, int numFuncs)
: grid_ (grid) , indexSet_(indexSet), numFuncs_ (numFuncs)
, myProc_ (grid.myProcessor() )
{
initialize();
}
//! send and recive DiscreteFunction
template <class DiscFuncType>
void sendRecive( DiscFuncType & vec)
{
#ifndef YGRID
loadSendBuffer(vec);
communicate(vec);
unloadRecvBuffer(vec);
#else
vec.print(std::cout,1);
typedef typename DiscFuncType::DofIteratorType DofIteratorType;
DofIteratorType it = vec.dbegin(-1);
GathScatt<double> fake;
grid_.communicate(it,fake,InteriorBorder_All_Interface,BackwardCommunication,grid_.maxlevel());
grid_.communicate(it,fake,Overlap_All_Interface,BackwardCommunication,grid_.maxlevel());
#endif
}
//! send and recive DiscreteFunction
template <class DiscFuncType>
void loadSendBuffer( DiscFuncType & vec, int pos=0)
{
#ifndef YGRID
{
typedef typename DiscFuncType::DofIteratorType DofIteratorType;
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
DofIteratorType it = vec.dbegin(-1);
for( ; procit != procend; ++procit)
{
(*procit).loadSendBuffer(it,pos);
}
}
#endif
}
//! send and recive DiscreteFunction
template <class DiscFuncType>
void unloadRecvBuffer( DiscFuncType & vec, int pos=0)
{
#ifndef YGRID
{
typedef typename DiscFuncType::DofIteratorType DofIteratorType;
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
DofIteratorType it = vec.dbegin(-1);
for( ; procit != procend; ++procit)
{
(*procit).unloadRecvBuffer(it);
}
}
#endif
}
//! send and recive DiscreteFunction
template <class DiscFuncList>
void sendReciveList( DiscFuncList & list)
{
assert(numFuncs_ == list.size());
{
typedef typename DiscFuncList::Iterator DLIterator;
DLIterator dl = list.begin();
DLIterator de = list.end();
int count = 0;
for(; dl != de ; ++dl )
{
loadSendBuffer(* (*dl) ,count );
count++;
}
}
communicate();
{
typedef typename DiscFuncList::Iterator DLIterator;
DLIterator dl = list.begin();
DLIterator de = list.end();
int count = 0;
for(; dl != de ; ++dl )
{
unloadRecvBuffer(* (*dl) , count);
count++;
}
}
}
// minimize timestepsize over all processors
template <typename T>
T timeStepSize(T timestep)
{
T ret=-1.0;
MPI_Allreduce(&timestep, &ret, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
return(ret);
}
// minimize timestepsize over all processors
void waitForAll()
{
MPI_Barrier ( MPI_COMM_WORLD );
}
// send and recive a vector
template <class DofArrayType>
void sendReciveVec( DofArrayType & vec)
{
#ifndef YGRID
{
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
for( ; procit != procend; ++procit)
{
(*procit).loadSendBuffer(vec);
}
}
communicate();
{
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
for( ; procit != procend; ++procit)
{
(*procit).unloadRecvBuffer(vec);
}
}
#endif
}
private:
// send and recive data to all neighbouring processors
void communicate()
{
#ifndef YGRID
// send data
{
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
//int tag = myProc_;
int no=0;
for( ; procit != procend; ++procit)
{
//std::cout << " Sende Proc " << (*procit).processor() << "\n";
MPI_Issend( (void *)(*procit).getSendBuffer(),(*procit).sendSize(),
MPI_DOUBLE, (*procit).processor(),
//tag , AlbertGridComm,
COMMUNICATOR_COMM_TAG, MPI_COMM_WORLD ,
(*procit).sendreq());
no++;
}
}
MPI_Status status;
// recive date
{
//int tag = myProc_;
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
int no=0;
for( ; procit != procend; ++procit)
{
MPI_Irecv( (void *)(*procit).getRecvBuffer(),(*procit).recvSize(),
MPI_DOUBLE, (*procit).processor(),
COMMUNICATOR_COMM_TAG, MPI_COMM_WORLD ,
(*procit).recvreq());
no++;
}
}
{
{
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
int no=0;
for( ; procit != procend; ++procit)
{
MPI_Wait((*procit).sendreq(),&status);
no++;
}
}
}
{
{
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
int no=0;
for( ; procit != procend; ++procit)
{
MPI_Wait((*procit).recvreq(),&status);
no++;
}
}
}
#endif
}
// build address book
void initialize ()
{
#ifndef YGRID
typedef typename GridType::template Traits<0>::LevelIterator LevelIteratorType;
LevelIteratorType it = grid_.template lbegin<0> (0,Ghosts);
LevelIteratorType endit = grid_.template lend <0> (0,Ghosts);
for( ; it != endit; ++it )
{
int owner = it->owner();
if( owner != myProc_ )
{
bool wehaveit = false;
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
for( ; procit != procend; ++procit)
{
if( (*procit).processor() == owner) wehaveit = true;
}
if(!wehaveit)
{
ProcListElement<double> tmp (owner, numFuncs_ );
otherProcs_.insert_before(otherProcs_.begin(), tmp);
}
}
}
// now we have the link list for the processors
typedef typename GridType::LeafIterator LeafIteratorType;
{
typedef typename ProcessorListType::Iterator ListIterator;
ListIterator procit = otherProcs_.begin();
ListIterator procend = otherProcs_.end();
for( ; procit != procend; ++procit)
{
{
LeafIteratorType it = grid_.leafbegin ( grid_.maxlevel(), Border, (*procit).processor() );
LeafIteratorType endit = grid_.leafend ( grid_.maxlevel(), Border, (*procit).processor() );
for(; it != endit; ++it)
{
(*procit).insertSend( indexSet_.template index<0> (*it, 0 ) );
}
}
{
LeafIteratorType it = grid_.leafbegin ( grid_.maxlevel(), Ghosts, (*procit).processor() );
LeafIteratorType endit = grid_.leafend ( grid_.maxlevel(), Ghosts, (*procit).processor() );
for(; it != endit; ++it)
{
(*procit).insertRecive( indexSet_.template index<0>(*it, 0 ) );
}
}
// malloc buffer
(*procit).makeBuffer();
//(*procit).print(std::cout);
}
}
#endif
}
// reference to corresponding grid
GridType & grid_;
// the index set
IndexSetType & indexSet_;
int numFuncs_;
// rank of my thread
const int myProc_;
// links to other processors
ProcessorListType otherProcs_;
};
#else
template <class GridType, class IndexSetType>
class AlbertGridCommunicator
{
public:
AlbertGridCommunicator(GridType &grid, IndexSetType &indexSet, int numFuncs) {}
template <class DiscFuncType>
void sendRecive( DiscFuncType & vec) {}
//! send and recive DiscreteFunction
template <class DiscFuncList>
void sendReciveList( DiscFuncList & list) {}
template <typename T>
T timeStepSize(T timestep)
{
return timestep;
}
};
#endif
template <class GridType, class CritType>
void makeParallelGrid (GridType &grid, CritType &crit)
{
for(int l=0; l <= grid.maxlevel(); l++)
{
typedef typename GridType::template Traits<0>::LevelIterator LevelIteratorType;
LevelIteratorType it = grid.template lbegin<0> (l);
LevelIteratorType endit = grid.template lend <0> (l);
for( ; it != endit; ++it )
{
crit.classify( *it );
}
}
}
} // end namespace Dune
#endif
grid/albertgrid/agelementindex.cc deleted 100644 → 0
+ 0
− 38
View file @ 4289ea0e
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef __ALBERTGRID_ELMEM_CC__
#define __ALBERTGRID_ELMEM_CC__
namespace AlbertHelp {
// IndexManagerType defined in albertgrid.hh
static IndexManagerType * tmpIndexStack = 0;
static void initIndexManager_elmem_cc(IndexManagerType * newIm)
{
tmpIndexStack = newIm;
assert(tmpIndexStack != 0);
}
static void removeIndexManager_elmem_cc()
{
tmpIndexStack = 0;
}
/* return the new element index for el->index */
static int get_elIndex()
{
assert(tmpIndexStack != 0);
return tmpIndexStack->getIndex();
}
/* when element is deleted remember the index */
static void free_elIndex(int ind)
{
assert(tmpIndexStack != 0);
tmpIndexStack->freeIndex(ind);
}
} // end namespace AlbertHelp
#endif
grid/albertgrid/agindex.hh deleted 100755 → 0
+ 0
− 259
View file @ 4289ea0e
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef __DUNE_AGINDEX_HH__
#define __DUNE_AGINDEX_HH__
#include "../../common/array.hh"
namespace Dune {
enum INDEXSTATE { NEW, OLD , USED, UNUSED };
template <class GridType>
class SerialIndexSet
{
typedef int T;
Array<int> globalIndex_;
Array<int> oldGlobalIndex_;
Array<INDEXSTATE> state_;
int maxIndex_;
int nextFreeIndex_;
int nextIndex_;
GridType & grid_;
public:
SerialIndexSet (GridType & grid) : grid_ (grid)
{
maxIndex_ = 0;
};
// calculate new highest index
//template <class GridType>
void insertNew (GridType & grid)
{
this->resize ( grid );
for(int l=0; l<=grid.maxlevel(); l++)
{
typedef typename GridType::template Traits<0>::LevelIterator LevelIteratorType;
LevelIteratorType endit = grid.template lend<0>(l);
LevelIteratorType it = grid.template lbegin<0>(l);
for( ; it != endit ; ++it)
{
this->insert( *it );
}
}
this->finish();
}
//template <class GridType>
void resize (GridType & grid)
{
this->resize ( grid.global_size (0));
}
// calculate new highest index
void resize (int newMaxInd )
{
if( globalIndex_.size() < newMaxInd )
{
globalIndex_.swap ( oldGlobalIndex_ );
globalIndex_.resize ( 2*newMaxInd );
state_.resize ( 2*newMaxInd );
for(int i=0; i<globalIndex_.size(); i++)
globalIndex_[i] = -1;
maxIndex_ = 2*newMaxInd;
nextFreeIndex_ = -1;
for(int i=0; i<oldGlobalIndex_.size(); i++)
{
globalIndex_[i] = oldGlobalIndex_[i];
if(state_[i] != UNUSED)
state_[i] = OLD;
if(globalIndex_[i] > nextFreeIndex_)
nextFreeIndex_ = globalIndex_[i];
if(globalIndex_[i] < 0)
state_[i] = UNUSED;
}
for(int i=oldGlobalIndex_.size(); i<state_.size(); i++)
state_[i] = UNUSED;
nextFreeIndex_++;
nextIndex_ = 0;
}
else
{
for(int i=0; i<state_.size(); i++)
{
if(state_[i] != UNUSED)
state_[i] = OLD;
}
}
//std::cout << nextFreeIndex_ << " freeInd \n";
}
void finish ()
{
for(int i=0; i<state_.size(); i++)
{
if(state_[i] == OLD)
{
state_[i] = UNUSED;
}
}
std::cout << maxIndex() << " max Index of Set \n";
}
int maxIndex () const
{
int max = 0;
for(int i =0; i<globalIndex_.size(); i++)
{
if(globalIndex_[i] > max ) max = globalIndex_[i];
}
return max;
}
int searchNext ()
{
if(nextIndex_ >= maxIndex_)
return -1;
while((state_[nextIndex_] != UNUSED) ||
(globalIndex_[nextIndex_] < 0) )
{
nextIndex_++;
if(nextIndex_ >= maxIndex_ )
return -1;
}
//std::cout << state_[nextIndex_] << " " << globalIndex_[nextIndex_] << "\n";
nextIndex_++;
return globalIndex_[nextIndex_-1];
}
// memorise index
template <class EntityType>
void insert (EntityType & en)
{
this->insert ( en.global_index() );
}
// memorise index
void insert (int num )
{
assert(num < globalIndex_.size() );
//std::cout << " insert num = " << num << "\n";
if(globalIndex_[num] < 0)
{
int ind = searchNext ();
if( ind >= 0 )
{
globalIndex_[num] = ind;
}
else
{
globalIndex_[num] = nextFreeIndex_;
nextFreeIndex_++;
}
state_[num] = NEW;
return;
}
state_[num] = USED;
}
void print ( ) const
{
std::cout << "Size " << globalIndex_.size() << "\n";
std::cout << "i | val | state \n";
for(int i=0; i<globalIndex_.size(); i++)
{
std::cout << i << " | " << globalIndex_[i] << " | " << state_[i] << "\n";
}
}
bool write_xdr(const char * filename, int timestep)
{
FILE *file;
XDR xdrs;
const char *path = NULL;
const char * fn = genFilename(path,filename, timestep);
file = fopen(fn, "wb");
if (!file)
{
printf( "\aERROR in AGIndexSet::write_xdr(..): couldnot open <%s>!\n", filename);
fflush(stderr);
return false;
}
xdrstdio_create(&xdrs, file, XDR_ENCODE);
this->processXdr(&xdrs);
xdr_destroy(&xdrs);
fclose(file);
}
bool read_xdr(const char * filename , int timestep)
{
FILE *file;
XDR xdrs;
const char *path = NULL;
const char * fn = genFilename(path,filename, timestep);
std::cout << "Reading <" << fn << "> \n";
file = fopen(fn, "rb");
if(!file)
{
printf( "\aERROR in AGIndexSet::read_xdr(..): couldnot open <%s>!\n", filename);
fflush(stderr);
return(false);
}
// read xdr
xdrstdio_create(&xdrs, file, XDR_DECODE);
this->processXdr(&xdrs);
xdr_destroy(&xdrs);
fclose(file);
return true;
}
bool processXdr(XDR *xdrs)
{
xdr_int ( xdrs, &maxIndex_ );
xdr_int ( xdrs, &nextIndex_ );
xdr_int ( xdrs, &nextFreeIndex_ );
globalIndex_.processXdr(xdrs);
state_.processXdr(xdrs);
return true;
}
int size () const
{
return nextFreeIndex_;
//return grid_.global_size(0);
}
bool isNew (int index) const
{
if(state_[index] == NEW)
return true;
return false;
}
int operator [] (int i) const
{
//printf(" gIndex_[%d] = %d \n",i,globalIndex_[i]);
return globalIndex_[i];
}
};
} // end namespace
#endif
grid/albertgrid/agmemory.hh deleted 100644 → 0
+ 0
− 173
View file @ 4289ea0e
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef __DUNE_AGMEMORY_HH__
#define __DUNE_AGMEMORY_HH__
namespace Dune
{
// organize the memory management for entitys used by the NeighborIterator
template <class Object>
class MemoryProvider
{
public:
typedef Object ObjectType;
struct ObjectEntity
{
ObjectEntity () : next (0), item (0) {};
ObjectEntity *next;
Object *item;
};
// freeEntity_ = NULL
MemoryProvider() : freeEntity_ (0) {};
// do not copy pointers
MemoryProvider(const MemoryProvider<Object> & org) : freeEntity_ (0) {}
// call deleteEntity
~MemoryProvider ();
// delete recursive all free ObjectEntitys
void deleteEntity(ObjectEntity *obj);
// i.e. return pointer to Entity
template <class GridType>
ObjectEntity *getNewObjectEntity(GridType &grid, int level);
// i.e. return pointer to Entity
template <class FuncSpaceType, class DofVecType>
ObjectEntity *getNewObjectEntity(const FuncSpaceType &f, DofVecType &d);
// i.e. get pointer to element
ObjectEntity * getNewObjectEntity();
// free, move element to stack, returns NULL
ObjectEntity * freeObjectEntity (ObjectEntity *obj);
private:
ObjectEntity *freeEntity_;
};
//************************************************************************
//
// MemoryProvider
//
//************************************************************************
template <class Object> template <class GridType>
inline typename MemoryProvider<Object>::ObjectEntity *
MemoryProvider<Object>::getNewObjectEntity
(GridType &grid, int level )
{
if(!freeEntity_)
{
ObjectEntity *oe = new ObjectEntity ();
oe->item = new Object (grid,level);
return oe;
}
else
{
ObjectEntity *oe = freeEntity_;
freeEntity_ = oe->next;
return oe;
}
}
template <class Object> template <class FuncSpaceType, class DofVecType>
inline typename MemoryProvider<Object>::ObjectEntity *
MemoryProvider<Object>::getNewObjectEntity(const FuncSpaceType &f , DofVecType &d )
{
if(!freeEntity_)
{
ObjectEntity *oe = new ObjectEntity ();
oe->item = new Object (f,d);
return oe;
}
else
{
ObjectEntity *oe = freeEntity_;
freeEntity_ = oe->next;
return oe;
}
}
template <class Object>
inline typename MemoryProvider<Object>::ObjectEntity *
MemoryProvider<Object>::getNewObjectEntity()
{
if(!freeEntity_)
{
ObjectEntity *oe = new ObjectEntity ();
oe->item = new Object ();
return oe;
}
else
{
ObjectEntity *oe = freeEntity_;
freeEntity_ = oe->next;
return oe;
}
}
template <>
inline MemoryProvider<ALBERT EL_INFO>::ObjectEntity *
MemoryProvider<ALBERT EL_INFO>::getNewObjectEntity()
{
if(!freeEntity_)
{
ObjectEntity *oe = new ObjectEntity ();
oe->item = (ALBERT EL_INFO *) std::malloc (sizeof(ALBERT EL_INFO));
return oe;
}
else
{
ObjectEntity *oe = freeEntity_;
freeEntity_ = oe->next;
return oe;
}
}
template <class Object>
inline MemoryProvider<Object>::~MemoryProvider()
{
if(freeEntity_) deleteEntity(freeEntity_);
}
template <class Object>
inline typename MemoryProvider<Object>::ObjectEntity *
MemoryProvider<Object>::freeObjectEntity(ObjectEntity *obj)
{
obj->next = freeEntity_;
freeEntity_ = obj;
return 0;
}
template <class Object>
inline void MemoryProvider<Object>::deleteEntity(ObjectEntity *obj)
{
if(obj)
{
if(obj->next)
deleteEntity(obj->next);
if(obj->item) delete obj->item;
delete obj;
}
}
template <>
inline void MemoryProvider<ALBERT EL_INFO>::deleteEntity(ObjectEntity *obj)
{
std::free(obj->item);
delete obj;
}
typedef MemoryProvider < ALBERT EL_INFO > ElInfoProvider;
static ElInfoProvider elinfoProvider;
} // end namespace Dune
#endif
grid/albertgrid/albertextra.hh deleted 100644 → 0
+ 0
− 1207
View file @ 4289ea0e
This diff is collapsed.
grid/albertgrid/albertgrid.cc deleted 100644 → 0
+ 0
− 4365
View file @ 4289ea0e
This diff is collapsed.
grid/albertgrid/part.cc deleted 100644 → 0
+ 0
− 490
View file @ 4289ea0e
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
/****************************************************************************/
/* ALBERT: an Adaptive multi Level finite element toolbox using */
/* Bisectioning refinement and Error control by Residual */
/* Techniques */
/* */
/* file: 2d/macro.c */
/* !!! includes Common/macro_common.c !!! */
/* */
/* description: tools for reading and writing macro triangulations */
/* */
/* */
/* history: */
/* */
/* */
/****************************************************************************/
/* */
/* authors: Alfred Schmidt and Kunibert G. Siebert */
/* Institut f"ur Angewandte Mathematik */
/* Albert-Ludwigs-Universit"at Freiburg */
/* Hermann-Herder-Str. 10 */
/* 79104 Freiburg */
/* Germany */
/* */
/* email: alfred,kunibert@mathematik.uni-freiburg.de */
/* */
/* (c) by A. Schmidt und K.G. Siebert (1996) */
/* */
/****************************************************************************/
/****************************************************************************/
/* modifications: vector oriented storage of macro triangulation data */
/* in structure MACRO_DATA; data is then converted into ALBERT-macro */
/* Author: Daniel Koester (2001) */
/****************************************************************************/
#include <ctype.h>
#include <rpc/types.h>
#include <rpc/xdr.h>
#include <algorithm>
#if LINKEN
#include <albert.h>
namespace Albert
{
#endif
/****************************************************************************/
/* these definitions are necessary at the moment for the macros "MEM_ALLOC" */
/* and "MEM_FREE" to work properly */
/****************************************************************************/
typedef int INT_NVERT[N_VERTICES];
typedef int INT_NNEIGH[N_NEIGH];
typedef S_CHAR SCHAR_NNEIGH[N_NEIGH];
typedef U_CHAR UCHAR_NNEIGH[N_NEIGH];
/****************************************************************************/
/* default_boundary(): This procedure returns a constant pointer to a */
/* BOUNDARY structure corresponding to a standard polygonal boundary. */
/* NOTE: The total number of possible boundary types is unknown during the */
/* first calls of this procedure. It therefore has to allocate a new */
/* BOUNDARY structure for each new argument "bound" passed to it; the */
/* address of this new BOUNDARY is stored in a dynamic array default_bounds.*/
/****************************************************************************/
const BOUNDARY *default_boundary(MESH *mesh, int bound)
{
FUNCNAME("default_boundary");
typedef BOUNDARY *BOUND_PTR;
static BOUND_PTR *default_bounds = nil;
static unsigned int i = 0, size = 1;
if(!bound) {
WARNING("tried to assign a BOUNDARY * with bound == 0!\n");
return(nil);
}
if(!default_bounds) {
default_bounds = MEM_ALLOC(size, BOUND_PTR);
default_bounds[0] = MEM_ALLOC(1, BOUNDARY);
default_bounds[0]->param_bound = nil;
default_bounds[0]->bound = bound;
return((const BOUNDARY *) default_bounds[0]);
}
for(i = 0; i < size; i++)
if (default_bounds[i]->bound == bound)
return((const BOUNDARY *) default_bounds[i]);
if(i == size) {
default_bounds = MEM_REALLOC(default_bounds, size, size + 1, BOUND_PTR);
default_bounds[size] = MEM_ALLOC(1, BOUNDARY);
default_bounds[size]->param_bound = nil;
default_bounds[size]->bound = bound;
size++;
}
return((const BOUNDARY *) default_bounds[size - 1]);
}
/****************************************************************************/
/* fill_more_boundary_dofs(): */
/* adds dof's at the edges of a given macro triangulation and calculates */
/* the number of edges */
/* function pointer bdry only used in 3d! */
/****************************************************************************/
static void fill_more_boundary_dofs(MESH *mesh,
const BOUNDARY *(*bdry)(MESH *, int))
{
int i, n_edges = 0, node = mesh->node[EDGE];
MACRO_EL *mel;
DOF *dof;
for (mel = mesh->first_macro_el; mel; mel = mel->next)
{
for (i = 0; i < N_NEIGH; i++)
{
if (!mel->neigh[i] || (mel->neigh[i]->index < mel->index))
{
n_edges++;
if (mesh->n_dof[EDGE])
{
dof = mel->el->dof[node+i] = get_dof(mesh, EDGE);
if (mel->neigh[i])
mel->neigh[i]->el->dof[node+mel->opp_vertex[i]] = dof;
}
}
}
}
mesh->n_edges = n_edges;
return;
}
/****************************************************************************/
/* fill_bound_info(): */
/* fills boundary information for the vertices of the macro triangulation */
/* The type of a boundary vertex is equal to the highest type of all */
/* adjoining boundary edges. If there are no boundary edges containing the */
/* vertex, it is assumed to be an interior vertex. */
/****************************************************************************/
static void fill_bound_info(MESH *mesh, MACRO_DATA *data,
const BOUNDARY *(*bdry)(MESH *, int))
{
FUNCNAME("fill_bound_info");
MACRO_EL *mel = mesh->first_macro_el;
int i, j, ne = mesh->n_elements, nv = mesh->n_vertices;
S_CHAR *bound = MEM_ALLOC(nv, S_CHAR);
for(i = 0; i < data->n_macro_elements; i++) {
for(j = 0; j < N_NEIGH; j++) {
if (data->boundary[i][j] != INTERIOR)
mel[i].boundary[j] = bdry(mesh, data->boundary[i][j]);
else
mel[i].boundary[j] = nil;
}
}
for (i = 0; i < nv; i++)
bound[i] = INTERIOR;
for (i = 0; i < ne; i++)
{
for (j = 0; j < N_NEIGH; j++)
{
if (mel[i].boundary[j])
{
if (mel[i].boundary[j]->bound >= DIRICHLET)
{
int j1 = data->mel_vertices[i][(j+1)%3];
int j2 = data->mel_vertices[i][(j+2)%3];
bound[j1] = std::max(bound[j1], mel[i].boundary[j]->bound);
bound[j2] = std::max(bound[j2], mel[i].boundary[j]->bound);
}
else if (mel[i].boundary[j]->bound <= NEUMANN)
{
int j1 = data->mel_vertices[i][(j+1)%3];
int j2 = data->mel_vertices[i][(j+2)%3];
if (bound[j1] != INTERIOR)
bound[j1] = std::max(bound[j1], mel[i].boundary[j]->bound);
else
bound[j1] = mel[i].boundary[j]->bound;
if (bound[j2] != INTERIOR)
bound[j2] = std::max(bound[j2], mel[i].boundary[j]->bound);
else
bound[j2] = mel[i].boundary[j]->bound;
}
}
}
}
for (i = 0; i < ne; i++)
for (j = 0; j < N_VERTICES; j++)
mel[i].bound[j] = bound[data->mel_vertices[i][j]];
MEM_FREE(bound, nv, S_CHAR);
return;
}
/****************************************************************************/
/* fill_best_edges(): The main job of this routine is to fill the arrays */
/* best_edges[] and neighs[] below. best_edges[elem] is best explained with */
/* some examples: */
/* best_edges[elem] == {2, 3, 3}: one longest edge, namely 2 */
/* best_edges[elem] == {0, 1, 3}: two longest edges, namely 0 and 1 */
/* best_edges[elem] == {2, 0, 1}: three longest edges, namely 2, 0, 1 */
/* neighs[elem] contains the global indices of the neighbour edges ordered */
/* by length to match best_edges[elem]. */
/****************************************************************************/
static void fill_best_edges(MACRO_DATA *data, int elem, UCHAR_NNEIGH edge,
INT_NNEIGH neighs)
{
static U_CHAR i;
static REAL l[3];
for(i = 0; i < N_EDGES; i++) {
l[i] = DIST_DOW(data->coords[data->mel_vertices[elem][(i + 1) % N_EDGES]],
data->coords[data->mel_vertices[elem][(i + 2) % N_EDGES]]);
edge[i] = i;
}
for (i = 0; i < N_EDGES; i++) {
if (l[i] > l[edge[0]]) edge[0] = i;
if (l[i] < l[edge[2]]) edge[2] = i;
}
edge[1] = N_EDGES - edge[0] - edge[2];
for(i = 0; i < N_NEIGH; i++)
neighs[i] = data->neigh[elem][edge[i]];
for (i = 1; i < N_EDGES; i++)
if ( (l[edge[i - 1]] - l[edge[i]]) > REAL_EPSILON * l[edge[i]] )
break;
for (; i < N_EDGES; i++)
edge[i]=N_EDGES;
}
/****************************************************************************/
/* new_refine_edge(): change the local indices of vertices, neighbours and */
/* boundaries to fit the choice of the new refinement edge. */
/* NOTE: the element's orientation is unimportant for this routine. */
/****************************************************************************/
static void new_refine_edge(MACRO_DATA *data, int elem, U_CHAR new_edge)
{
static int buffer_vertex, buffer_neigh;
static S_CHAR buffer_boundary;
if(new_edge != 2) {
buffer_vertex = data->mel_vertices[elem][0];
buffer_neigh = data->neigh[elem][0];
buffer_boundary = data->boundary[elem][0];
switch (new_edge) {
case 0 :
data->mel_vertices[elem][0] = data->mel_vertices[elem][1];
data->mel_vertices[elem][1] = data->mel_vertices[elem][2];
data->mel_vertices[elem][2] = buffer_vertex;
data->neigh[elem][0] = data->neigh[elem][1];
data->neigh[elem][1] = data->neigh[elem][2];
data->neigh[elem][2] = buffer_neigh;
data->boundary[elem][0] = data->boundary[elem][1];
data->boundary[elem][1] = data->boundary[elem][2];
data->boundary[elem][2] = buffer_boundary;
break;
case 1 :
data->mel_vertices[elem][0] = data->mel_vertices[elem][2];
data->mel_vertices[elem][2] = data->mel_vertices[elem][1];
data->mel_vertices[elem][1] = buffer_vertex;
data->neigh[elem][0] = data->neigh[elem][2];
data->neigh[elem][2] = data->neigh[elem][1];
data->neigh[elem][1] = buffer_neigh;
data->boundary[elem][0] = data->boundary[elem][2];
data->boundary[elem][2] = data->boundary[elem][1];
data->boundary[elem][1] = buffer_boundary;
}
}
}
/****************************************************************************/
/* reorder(): the main routine for correcting cyles */
/* */
/* Refinement edges are chosen using the following priority: */
/* 0. If the current element elem only has one longest edge then */
/* choose that edge as refinement edge. */
/* 1. If possible, choose the refinement edge as one of the longest */
/* edges along the boundary. */
/* 2. Otherwise chose the refinement edge as one of the longest edges */
/* whose corresponding neighbour's edge is also one of its longest */
/* ones (thus creating compatibly divisible pairs of elements) */
/* 3. Choose a longest edge towards an already tested element. */
/* 4. If all else fails, choose an arbitrary longest edge. */
/****************************************************************************/
static void reorder(MACRO_DATA *data, U_CHAR *test, int elem,
INT_NNEIGH *neighs, UCHAR_NNEIGH *best_edges)
{
FUNCNAME("reorder");
static U_CHAR j, k;
MSG("Current elem: %d, best_edges: %d %d %d\n", elem, best_edges[elem][0], best_edges[elem][1], best_edges[elem][2]);
test[elem] = true;
if (best_edges[elem][1] == N_EDGES) {
new_refine_edge(data, elem, best_edges[elem][0]);
return;
}
for (j = 0; best_edges[elem][j] < N_EDGES; j++) {
MSG("Looking at best_edges[%d][%d]...\n", elem, j);
if (neighs[elem][j] < 0) {
MSG("It is a border edge! Selecting it...\n");
new_refine_edge(data, elem, best_edges[elem][j]);
return;
}
if (!test[neighs[elem][j]]) {
for (k = 0; best_edges[neighs[elem][j]][k] < N_EDGES; k++)
if(neighs[neighs[elem][j]][k] == elem) {
MSG("Found compatibly divisible neighbour %d!\n", neighs[elem][j]);
test[neighs[elem][j]] = true;
new_refine_edge(data, elem, best_edges[elem][j]);
new_refine_edge(data,neighs[elem][j],best_edges[neighs[elem][j]][k]);
return;
}
}
}
MSG("No immediate patch found - trying to select an edge towards tested elements.\n");
for (j = 0; best_edges[elem][j] < N_EDGES; j++) {
MSG("Looking at best_edges[%d][%d]...\n", elem, j);
if (test[neighs[elem][j]]) {
MSG("Found tested neighbour on edge %d.", j);
new_refine_edge(data, elem, best_edges[elem][j]);
return;
}
}
MSG("Finally resorted to selecting edge %d.\n", best_edges[elem][0]);
new_refine_edge(data, elem, best_edges[elem][0]);
return;
}
/****************************************************************************/
/* correct_cycles(): Correct refinement cycles using reorder() */
/****************************************************************************/
static void correct_cycles(MACRO_DATA *data)
{
FUNCNAME("correct_cycles");
int elem;
U_CHAR *test;
INT_NNEIGH *neighs = MEM_ALLOC(data->n_macro_elements, INT_NNEIGH);
UCHAR_NNEIGH *best_edges = MEM_ALLOC(data->n_macro_elements, UCHAR_NNEIGH);
test = MEM_CALLOC(data->n_macro_elements, U_CHAR);
for(elem = 0; elem < data->n_macro_elements; elem++)
fill_best_edges(data, elem, best_edges[elem], neighs[elem]);
for(elem = 0; elem < data->n_macro_elements; elem++)
if(!test[elem]) reorder(data,test,elem, neighs, best_edges);
MEM_FREE(test, data->n_macro_elements, U_CHAR);
MEM_FREE(neighs, data->n_macro_elements, INT_NNEIGH);
MEM_FREE(best_edges, data->n_macro_elements, UCHAR_NNEIGH);
}
/****************************************************************************/
/* orientation(): checks and corrects whether the element is oriented */
/* counterclockwise. */
/****************************************************************************/
U_CHAR orientation(MACRO_DATA *data)
{
int i, vert_buffer, neigh_buffer;
REAL_D e1, e2;
REAL det, *a0;
S_CHAR bound_buffer;
U_CHAR result = false;
for(i = 0; i < data->n_macro_elements; i++) {
a0 = data->coords[data->mel_vertices[i][0]];
e1[0] = data->coords[data->mel_vertices[i][1]][0] - a0[0];
e1[1] = data->coords[data->mel_vertices[i][1]][1] - a0[1];
e2[0] = data->coords[data->mel_vertices[i][2]][0] - a0[0];
e2[1] = data->coords[data->mel_vertices[i][2]][1] - a0[1];
det = e1[0]*e2[1] - e1[1]*e2[0];
if(det < 0) {
result = true;
vert_buffer = data->mel_vertices[i][0];
data->mel_vertices[i][0] = data->mel_vertices[i][1];
data->mel_vertices[i][1] = vert_buffer;
neigh_buffer = data->neigh[i][0];
data->neigh[i][0] = data->neigh[i][1];
data->neigh[i][1] = neigh_buffer;
bound_buffer = data->boundary[i][0];
data->boundary[i][0] = data->boundary[i][1];
data->boundary[i][1] = bound_buffer;
}
}
return(result);
}
/****************************************************************************/
/* macro_test(): Check data for potential cycles during refinement */
/* At the moment, a correction (and subsequent writing of a correct macro */
/* data file) can only be done in 2D */
/* Additionally, the element orientation is checked and corrected. */
/****************************************************************************/
static int cycles(MACRO_DATA *);
extern void macro_test(MACRO_DATA *data, const char *nameneu)
{
FUNCNAME("macro_test");
U_CHAR error_found = false;
int i = -1;
i = cycles(data);
if (i >= 0)
{
error_found = true;
WARNING("There is a cycle beginning in macro element %d.\n", i);
MSG("Correcting refinement edges....\n");
correct_cycles(data);
}
#if DIM_OF_WORLD == 2
if(orientation(data)) {
error_found = true;
WARNING("Element orientation was corrected for some elements.\n");
}
#endif
if (error_found && nameneu) {
MSG("Attempting to write corrected macro data to file %s...\n", nameneu);
//write_macro_data(data, nameneu);
}
return;
}
#include "partialGrid.cc"
#if LINKEN
} // namespace Albert
#endif
grid/albertgrid/partialGrid.cc deleted 100644 → 0
+ 0
− 1772
View file @ 4289ea0e
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment