FEM.hpp

 
//
// $Id: hpp,v 1.0 2000/01/01 00:00:01 werner Exp $
//
// $Log: hpp,v $
// Created 2012/05/16 16:22:16  marcel
// Initial version
//
//
#ifndef __FIBER_GRID_TYPES_FEM_HPP
#define __FIBER_GRID_TYPES_FEM_HPP
 
#include "gridtypesDllApi.h"
#include <grid/Grid.hpp>
#include <eagle/PhysicalSpace.hpp>
 
#include <aerie/KDTree.hpp>
#include <aerie/GeometricFunctions.hpp>
#include <aerie/BoundingBox.hpp>
 
#include <field/OnDemandCreator.hpp>
//#include "FEMFields.hpp"
 
 
//#define VERBOSE
 
namespace Fiber
{
 
typedef uint32_t FEMIndex_t;
typedef std::vector<FEMIndex_t> FEMIndices_t;
 
typedef std::vector<Eagle::PhysicalSpace::tvector> FEMDirections_t;
 
typedef std::vector<std::vector<FEMIndex_t> > FEMFaceList_t;
 
 
template<>
inline  std::string     element_to_string(const FEMFaceList_t&EdgeList)
{
string  retval = "{(" + std::to_string(EdgeList.size()) + ") ";
 
        if (EdgeList.size()>0)
                retval += "<>"; // std::to_string(EdgeList[0]);
 
        for(size_t i = 1; i<EdgeList.size(); i++)
        {
                retval += ", <>"; // + std::to_string(EdgeList[i]);
        }
 
        return retval + " }";
}
 
template<>
inline  std::string     element_to_string(const FEMDirections_t&FEM)
{
string  retval = "{(" + std::to_string(FEM.size()) + ") ";
 
        if (FEM.size()>0)
                retval += "<>"; // std::to_string(EdgeList[0]);
 
        for(size_t i = 1; i<FEM.size(); i++)
        {
                retval += ", <>"; // + std::to_string(EdgeList[i]);
        }
 
        return retval + " }";
}
 
 
 
enum {FEM_VERTS=0, FEM_CELLS};
 
namespace FEMFunc
{
typedef Eagle::PhysicalSpace::point   point_t;
typedef Eagle::PhysicalSpace::tvector tvector_t;
 
extern gridtypes_API
point_t getLocalC3D8GaussCoords( size_t IPNr );
extern gridtypes_API
point_t getLocalC3D20GaussCoords( size_t IPNr );
#ifdef NOT_IMPLEMENTED
point_t getLocalC3D20RGaussCoords( unsigned IPNr );
#endif
 
extern gridtypes_API
double getC3D8Weight( const point_t& local, size_t node);
extern gridtypes_API
double getC3D20Weight( const point_t& local, size_t node);
 
extern gridtypes_API
double getC3D8GaussWeight(   size_t gauss_nr, size_t node_nr );
extern gridtypes_API
double getC3D20GaussWeight(  size_t gauss_nr, size_t node_nr );
extern gridtypes_API
double getC3D20RGaussWeight( size_t gauss_nr, size_t node_nr );
 
 
 
 
// maybe better use a const array instead modulos for the -1 1 0 shuffling
// for vlen
template<class T>
T getInterpolatedC3D8( const point_t& local, const FEMIndices_t& coords_cell, const MultiArray<1, T>& data )
{
T ret;
        for( unsigned i = 0; i < T::Dims; i++ )
                ret[i] = 0;
 
 
        for( unsigned j = 0; j < 8; j++ )
        {
        double weight = getC3D8Weight( local, j );
 
                for( unsigned i = 0; i < T::Dims; i++ )
                        ret[i] += weight * T(data[ coords_cell[j] ])[i];
        }
        return ret;
}
 
 
 
//for fixed size
template<class T>
T getInterpolatedC3D8( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t l, const MultiArray<1, T>& data )
{
T ret;
        for( unsigned i = 0; i < T::Dims; i++ )
                ret[i] = 0;
 
MultiIndex<2> mi;
 
        for( unsigned j = 0; j < 8; j++ )
        {
                mi = j, l;
        double weight = getC3D8Weight( local, j );
 
                for( unsigned i = 0; i < T::Dims; i++ )
                        ret[i] += weight * T(data[ coords_cell[mi] ])[i];
        }
        return ret;
}
 
// specialization for non compound types
template<> gridtypes_API
float getInterpolatedC3D8<float>( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t l, const MultiArray<1, float>& data );
 
template<> gridtypes_API
double getInterpolatedC3D8<double>( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t l, const MultiArray<1, double>& data );
 
 
template<class T, int N>
T getInterpolatedC3D8Dn( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t cell, const MultiArray<1, T>& data )
{
T ret;
 Verbose(0) << "getInterpolatedC3D8Dn() derivation in N=" << N << " not implemented!";
 return ret;
}
 
template<> gridtypes_API
point_t getInterpolatedC3D8Dn<point_t, 0>( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t cell, const MultiArray<1, point_t>& data );
 
template<> gridtypes_API
point_t getInterpolatedC3D8Dn<point_t, 1>( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t cell, const MultiArray<1, point_t>& data );
 
template<> gridtypes_API
point_t getInterpolatedC3D8Dn<point_t, 2>( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t cell, const MultiArray<1, point_t>& data );
 
 
// maybe better use a const array instead modulos for the -1 1 0 shuffling
// for vlen
template<class T>
T getInterpolatedC3D20( const point_t& local, const FEMIndices_t& coords_cell, const MultiArray<1, T>& data )
{
T ret;
        for( unsigned i = 0; i < T::Dims; i++ )
                ret[i] = 0;
 
double weight;
        // corner nodes
        for(size_t j = 0; j < 20; j++)
        {
                weight = getC3D20Weight( local, j );
 
                for( unsigned i = 0; i < T::Dims; i++ )
                        ret[i] += weight * T(data[ coords_cell[j] ])[i];
        }
 
        return ret;
}
 
 
 
//for fixed size
template<class T>
T getInterpolatedC3D20( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t l, const MultiArray<1, T>& data )
{
T ret;
        for( unsigned i = 0; i < T::Dims; i++ )
                ret[i] = 0;
 
MultiIndex<2> mi;
 
double weight;
        // corner nodes
        for(size_t j = 0; j < 20; j++)
        {
                mi = j, l;
                weight = getC3D20Weight( local, j );
 
                for( unsigned i = 0; i < T::Dims; i++ )
                        ret[i] += weight * T(data[ coords_cell[mi] ])[i];
        }
 
        return ret;
}
 
// specialization for non compound types
template<>
float getInterpolatedC3D20<float>( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t l, const MultiArray<1, float>& data );
 
template<>
double getInterpolatedC3D20<double>( const point_t& local, const MultiArray<2, FEMIndex_t>& coords_cell, index_t l, const MultiArray<1, double>& data );
 
}
 
 
 
 
 
//namespace fem
//{
 
 
 
class gridtypes_API FEM
{
        enum { UNKNOWN, MIXED, C3D8, C3D8R, C3D20, C3D20R, C3D4, C3D10, C3D6, C3D15, T3D2,
               S2D4 } fe_type;
 
FEMFaceList_t face_list;
 
bool valid;
 
// use 32bit for opengl
RefPtr<Grid> v_grid;
RefPtr<Grid> i_grid;
 
RefPtr<Skeleton> v_vertices, v_cells;
RefPtr<Skeleton> i_vertices, i_cells;
 
RefPtr<Representation> v_cells_as_vertices, i_cells_as_vertices,
                       v_points, i_points, v_vertices_as_cells;
 
//must have
int cells_dim;
 
RefPtr<Field> v_coords;
RefPtr<Eagle::KDTree<3, index_t> > v_coords_tree;
 
RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > v_coords_arr; // <-- really load data here?
RefPtr<Field> v_cell_positions;
RefPtr<	MemArray<2, FEMIndex_t> > v_cell_arr;                   // <-- really load data here?
 
//RefPtr<Field> v_displacement;
RefPtr<Field> v_as_cell_positions;
 
 
RefPtr<Field> i_cell_positions;
//RefPtr<Field> i_distorsion;
MultiIndex<1> i_size;
//RefPtr<Field> i_tension;
//RefPtr<MemArray<1, Eagle::metric33> > i_tension_arr;
 
//optional (flags call of setup standard field (if null));
RefPtr<Field> i_coords;
 
public:
 
RefPtr<Field> getVertexField( string field_name )
{
        if(!v_grid)
                return NullPtr();
 
        return (*v_grid)(field_name);
}
 
RefPtr<Field> getDisplacementField( )
{
RefPtr<Field> displacement = getVertexField( "U" );
 
        if(!displacement)
        {
        FindPossibleDisplacementField displacement_finder( v_points );
                v_points->iterate( displacement_finder );
 
                displacement = displacement_finder.getBestField();
        }
 
        return displacement;
}
 
 
RefPtr<Field> getVertexCellField( string field_name )
{
        if(!v_cells_as_vertices)
                return NullPtr();
 
        return (*v_cells_as_vertices)(field_name);
}
 
 
RefPtr<Field> getIntegrationPointField( string field_name )
{
        if(!i_grid)
                return NullPtr();
 
        return (*i_grid)(field_name);
}
 
 
static  SkeletonID cellID3D()
{
return SkeletonID(3,2);
}
 
static  SkeletonID cellID2D()
{
return SkeletonID(2,2);
}
 
 
 
static double maxCellSize( RefPtr<MemArray<1, FEMIndices_t> >& cells );
 
template<class PointType>
static  RefPtr<Eagle::KDTree<PointType::SIZE, index_t> >  getCellMidPointTree( RefPtr<MemArray<1, FEMIndices_t> >& cell_set, RefPtr<MemArray<1, PointType> >& VertexSet );
 
RefPtr<MemArray<1, FEMIndices_t> > getIntPointCellPositions( const RefPtr<FragmentID>& f = NullPtr()  )
{
        if(!i_cell_positions)
                return NullPtr();
        return i_cell_positions->getData( f );
}
 
RefPtr<MemArray<1, FEMIndices_t> > getVertexCellPositions( const RefPtr<FragmentID>& f = NullPtr() )
{
        if(!v_cell_positions)
                return NullPtr();
 
        return v_cell_positions->getData( f );
}
 
MultiIndex<1> getIntPointSize()
{
        return  i_size;
}
 
MultiIndex<1> getVertexSize()
{
        if( !v_coords )
                return 0;
 
//Having no size interface does NOT mean there are no vertices!
        if (RefPtr<SizeInterface> SI = getSize( v_coords->getCreator() ))
        {
                return  SI->nElements();
        }
 
        return 0;
}
 
//@todo interface? not load cell array in constructor?
MultiIndex<2> getVertexCellSize()
{
        if( !v_cell_arr )
                return MultiIndex<2>(0,0);
 
        return  v_cell_arr->Size();
}
 
 
MultiIndex<2> getCellsSize()
{
        if( !v_cell_arr )
                return MultiIndex<2>(0,0);
 
        return v_cell_arr->Size();
}
// no fragment support here ... yet
RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > getVertexPositions()
{
        assert( v_coords_arr );
        return v_coords_arr;
}
 
std::pair<Eagle::PhysicalSpace::tvector, Eagle::PhysicalSpace::tvector> getCellExtendsInBB( const Eagle::BoundingBox& bb );
 
bool isValid()
{
        return valid;
}
 
struct  gridtypes_API PredefinedFieldNames
        {
        static  const char
                DISPLACEMENT[],
                MEAN_STRESS[],
                DEVIATORIC_RADIUS[],
                STRESS_INVARIANT1[],
 
                STRESS_INVARIANT2[],
 
                STRESS_INVARIANT3[],
                STRESS_MISES[],
                DEVIATORIC_STRESS[],
                STRESS[],
                AVERAGED[],
                VERT_AS_CELLS[],
                POSTFIX_IP[],
                POSTFIX_TRIANGULAR[];
};
 
        PredefinedFieldNames    FieldNames;
 
 
        FEM( const WeakPtr<Grid>&VG, const WeakPtr<Grid>&IG );
 
        ~FEM();
 
        struct gridtypes_API GetFirstFieldSize : public FieldIterator
        {
                MultiIndex<1> m_size;
 
                GetFirstFieldSize();
                bool apply (const FieldID &d, const Field &R);
        };
 
        struct gridtypes_API FindPossibleStressField : public FieldIterator
        {
                RefPtr<Representation> m_i_points;
                RefPtr<Field> m_i_tension;
                RefPtr<Field> m_first;
 
                FindPossibleStressField( const RefPtr<Representation> i_points );
                bool apply (const FieldID &d, const Field &R);
                RefPtr<Field> getBestField();
        };
 
        struct gridtypes_API FindPossibleDisplacementField : public FieldIterator
        {
                RefPtr<Representation> m_v_points;
                RefPtr<Field> m_v_displacement;
                RefPtr<Field> m_first;
 
                FindPossibleDisplacementField( const RefPtr<Representation> v_points );
                bool apply (const FieldID &d, const Field &R);
                RefPtr<Field> getBestField();
        };
 
        struct gridtypes_API SetupDerivedFields : public FieldIterator
        {
                RefPtr<Grid> m_v_grid, m_i_grid;
                RefPtr<Representation> m_i_points;
                MultiIndex<1> m_size;
 
                SetupDerivedFields( const RefPtr<Grid> v_grid, const RefPtr<Grid> i_grid, const RefPtr<Representation> i_points, MultiIndex<1> size );
                bool apply (const FieldID &d, const Field &R);
 
                template<class T>
                void setFieldCreator( string source_field, string field_name );
        };
 
         struct gridtypes_API ComputeToVertexCellFields : public FieldIterator
        {
                RefPtr<Grid> m_v_grid, m_i_grid;
                RefPtr<Representation> m_v_rep;
                MultiIndex<2> m_size;
 
                ComputeToVertexCellFields( const RefPtr<Grid> v_grid, const RefPtr<Grid> i_grid, const RefPtr<Representation> v_points, MultiIndex<2> size );
                bool apply (const FieldID &d, const Field &R);
 
                template<class T>
                void setFieldCreator( string source_field, string field_name );
        };
 
 
         struct gridtypes_API AverageToVertexFields : public FieldIterator
        {
                RefPtr<Grid> m_v_grid, m_i_grid;
                RefPtr<Representation> m_v_points;
                MultiIndex<1> m_size;
 
                AverageToVertexFields( const RefPtr<Grid> v_grid, const RefPtr<Grid> i_grid, const RefPtr<Representation> v_points, MultiIndex<1> size );
                bool apply (const FieldID &d, const Field &R);
 
                template<class T>
                void setFieldCreator( string source_field, string field_name );
        };
 
 
        RefPtr<Field> setUpStandardFields();
 
        bool isC3D8()  const;
        bool isC3D8R() const;
        bool isC3D20() const;
        bool isC3D20R()const;
        bool isMixed() const ;
        bool isUnknown() const;
        bool isC3D4()  const;
        bool isC3D10() const;
        bool isC3D6()  const;
        bool isC3D15() const;
        bool isT3D2()  const;
        bool isS2D4() const;
        // ...
 
        bool isHexa()  const;
        bool isTetra() const;
        bool isWedge() const;
        bool isBar()   const;
 
        RefPtr<Field> getVertexCells()
        {
                return v_cell_positions;
        }
 
 
        FEMFaceList_t& getFaceList()
        {
                return face_list;
        }
 
//      static RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > getLocalC3D8IPCoords();
 
 
 
 
/*
        RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > getVertexCoords( const RefPtr<FragmentID>& f = NullPtr() )
        {
                if( !v_coords )
                {
                        puts("FEm::getVertexCoords no vertex coords");
                        return NullPtr();
                }
                return v_coords->getData( f );
        }
 
        RefPtr<MemArray<1, FEMIndices_t> > getCellsAsVertices( const RefPtr<FragmentID>& f = NullPtr() )
        {
                if(!v_cell_positions)
                {
                        puts("FEm::getCellsAsVertices no vertex cells");
                        return NullPtr();
                }
                return v_cell_positions->getData( f );
        }
 
 
        RefPtr<MemArray<1, FEMIndices_t> > getCellsAsIntPoints( const RefPtr<FragmentID>& f = NullPtr() )
        {
                if(!i_cell_positions)
                {
                        puts("FEm::getCellsAsVertices no vertex cells");
                        return NullPtr();
                }
                return i_cell_positions->getData( f );
        }
*/
 
void getCellCandidates( std::set<index_t>&res_candidates,
                        const Eagle::PhysicalSpace::tvector& range, const Eagle::PhysicalSpace::point& center )
{
        res_candidates.clear();
 
//MultiArray< 2, FEMIndex_t >&   cells = *v_cell_arr;
//MultiArray<1, Eagle::point3>& crds = *v_coords_arr;
 
        if( !v_coords_tree )
        {
                Verbose(0) << "FEM::getCellCandidates() no tree!";
                exit(0);
        }
 
        if( ! v_as_cell_positions )
        {
                Verbose(0) << "FEM::getCellCandidates() vertices as cell array!";
                exit(0);
        }
 
RefPtr<MemArray<1, FEMIndices_t> > v_as_cell_arr;
        v_as_cell_arr = v_as_cell_positions->getData();
 
        assert( v_as_cell_arr );
 
MultiArray<1, FEMIndices_t>& v_cells = *v_as_cell_arr;
 
std::multimap<double, index_t> result_vertices;
 
        v_coords_tree->getManhattenRange( center, range, result_vertices );
 
//      assert( result_vertices.size() > 0 );
 
#ifdef VERBOSE
        std::cout << "found " << result_vertices.size() << " vertices via tree" << std::endl;
#endif
 
 
 
std::multimap<double, index_t>::const_iterator res_it;
        for( res_it = result_vertices.begin(); res_it != result_vertices.end(); res_it++ )
        {
        FEMIndices_t::const_iterator cell_it;
                for( cell_it = v_cells[res_it->second].begin() ; cell_it != v_cells[res_it->second].end(); cell_it++ )
                        res_candidates.insert( *cell_it );
//              break;
        }
 
#ifdef VERBOSE
std::set<index_t>::const_iterator cand_it;
        std::cout << "cell canditates are: " ;
        for( cand_it = res_candidates.begin(); cand_it != res_candidates.end(); cand_it++ )
                std::cout << *cand_it << " ";
        std::cout << std::endl;
#endif
}
 
void getCellCandidates( std::set<index_t>&res_candidates, const Eagle::BoundingBox&bb )
{
Eagle::PhysicalSpace::tvector range;
 
        range = ( bb.maxcoord() - bb.mincoord() ) / 2;
 
        getCellCandidates( res_candidates, range, bb.getCenter() );
}
 
 
 
// think about templatification
// prototype does it projectively
int getLocalCoordinate(  Eagle::point3& local, FEMIndex_t& cell, const Eagle::point3& world )
{
//RefPtr<       MemArray<2, FEMIndex_t> > cell_arr = v_cell_positions->getData();
//      assert(cell_arr);
 
std::set<index_t> cell_candidates;
 
Eagle::PhysicalSpace::tvector range;
        range.x() = 2;
        range.y() = 2;
        range.z() = 2;
 
 
       // todo: replacve by adjusted number
 
       getCellCandidates(cell_candidates, range, world );
 
MultiArray<1, Eagle::point3>& crds = *v_coords_arr;
MultiArray< 2, FEMIndex_t >&   cells = *v_cell_arr;
 
std::set<index_t>::const_iterator cand_it;
        for( cand_it = cell_candidates.begin(); cand_it != cell_candidates.end(); cand_it++ )
        {
 
 
                // project in candidates
        MultiIndex<2> max = v_cell_arr->Size();
 
                if( max[0] == 8 || max[0] == 20)
                        cuboidToUnitCube( local, world,
                                          crds[ cells[*cand_it][0] ],
                                          crds[ cells[*cand_it][4] ],
                                          crds[ cells[*cand_it][3] ],
                                          crds[ cells[*cand_it][1] ],
                                          crds[ cells[*cand_it][6] ] );
 
#ifdef VERBOSE
                std::cout << "getLocalCoordinate() checking cand: " << *cand_it << " : " << local << std::endl;
#endif
 
                if( local[0] >= 0.0 && local[1] >= 0.0 && local[2] >= 0.0 &&
                    local[0] <= 1.0 && local[1] <= 1.0 && local[2] <= 1.0 )
                {
#ifdef VERBOSE
                        std::cout << "point found in cell: " << *cand_it << std::endl;
#endif
                        // compute local crds from 0..1 to -1..+1
                        local += Eagle::PhysicalSpace::tvector(local) - Eagle::PhysicalSpace::tvector(1.,1.,1.);
                        //local[0] *= -1;
                        //local[1] *= -1;
                        //local[2] *= -1;
                        cell = FEMIndex_t(*cand_it);
 
#ifdef VERBOSE
                        std::cout << "getLocalCoordinate() EXIT 0" << std::endl;
                        std::cout.flush();
#endif
 
                        return 0;
                }
        }
 
 
#ifdef VERBOSE
        std::cout << "getLocalCoordinate() EXIT 1 " << std::endl;
        std::cout.flush();
#endif
        return 1;
}
 
//(  Eagle::point3& local, FEMIndex_t& cell, const Eagle::point3& world )
bool getLocalCoordinateNewton( Eagle::PhysicalSpace::point & local, const FEMIndex_t& cell,
                               const Eagle::PhysicalSpace::point&world, const double prec = 0.01, const unsigned max_iterations = 20 )
{
//      Verbose(0) << "getLocalCoordinateNewton() cell: " << cell;
//      Verbose(0) << "getLocalCoordinateNewton() nodes size: " << v_coords_arr->Size();
 
        if( cell > v_coords_arr->Size()[0] )
        {
                Verbose(0) << "getLocalCoordinateNewton() cell out of bounds!";
                return 1;
//              exit(0);
        }
 
typedef  Eagle::PhysicalSpace::point point;
typedef  Eagle::PhysicalSpace::tvector tvector;
 
unsigned iterations = 0;
point   local_iterative, world_iterative;
tvector dPdU, dPdV, dPdW;
tvector delta_world;
 
Eagle::Quadratic<3, double> sys, sys_inv;
Eagle::Column<3, double>    d_local, b;
 
// guess in mid of cell
        local_iterative = { 0.0, 0.0, 0.0 };
 
        world_iterative =  FEMFunc::getInterpolatedC3D8<point>( local_iterative, *v_cell_arr, cell, *v_coords_arr );
 
        // do newton iteration until a certain precision
        while( norm2(world - world_iterative) > prec*prec )
        {
                Verbose(0) << "iteration: " << iterations;
 
          dPdU = tvector( FEMFunc::getInterpolatedC3D8Dn<point, 0> ( local_iterative, *v_cell_arr, cell, *v_coords_arr ) );
          dPdV = tvector( FEMFunc::getInterpolatedC3D8Dn<point, 1> ( local_iterative, *v_cell_arr, cell, *v_coords_arr ) );
          dPdW = tvector( FEMFunc::getInterpolatedC3D8Dn<point, 2> ( local_iterative, *v_cell_arr, cell, *v_coords_arr ) );
 
          Verbose(60) << "dPdU: " << dPdU;
          Verbose(60) << "dPdV: " << dPdV;
          Verbose(60) << "dPdW: " << dPdW;
 
// solve system for (dU, dV, dW):
                // the 3x3 system:
                //    dPdU[0]*dU + dPdV[0]*dV + dPdW[0]*dW = delta_world[0]
                //    dPdU[1]*dU + dPdV[1]*dV + dPdW[1]*dW = delta_world[1]
                //    dPdU[2]*dU + dPdV[2]*dV + dPdW[2]*dW = delta_world[2]
 
                sys = dPdU[0], dPdU[1], dPdU[2],
                      dPdV[0], dPdV[1], dPdV[2],
                      dPdW[0], dPdW[1], dPdW[2];
 
                delta_world = world - world_iterative;
 
                Verbose(0) << "delta world: " << delta_world;
 
                try
                {
                        ComputeInverse(sys_inv, sys);
                }
                catch(const Eagle::DegeneratedMatrix&)
                {
                        Verbose(0) << "getLocalCoordinateNewton() Degenerated matrix error in inversion!";
                        assert(0);
                }
 
        Eagle::Matrix<3,3, double> &sysm_inv = sys_inv ;
        Eagle::Matrix<3,3, double> &sysm = sys;
 
                b = delta_world;
                d_local = sysm_inv * b;
 
        tvector d_world =  tvector::DomainVector_t(sysm * d_local);
 
                local_iterative[0] += d_local[0];
                local_iterative[1] += d_local[1];
                local_iterative[2] += d_local[2];
 
                world_iterative += d_world;
 
                Verbose(51) << "getLocalCoordinateNewton() in loop[" << iterations << "]: local=" << local_iterative;
 
                Verbose(0) <<"world: " << world;
                Verbose(0) <<"world_iterative: " << world_iterative;
                Verbose(0) <<"residuum: " << norm2(world - world_iterative);
//              local = local_iterative;
 
 
                if( iterations++ > max_iterations )
                {
                        Verbose(0) << "getLocalCoordinateNewton prevented endless loop";
                        return 1;
                }
        }
 
        local = local_iterative;
 
 
        for( unsigned i = 0; i < 3; i++ )
        {
                if( local_iterative[i] < -1 )
                {
                        Verbose(0) << "getLocalCoordinateNewton local out of -bounds: " << local_iterative;
//                              local_iterative[i] = -0.99;
                        return 1;
                }
                if( local_iterative[i] > 1. )
                {
                        Verbose(0) << "getLocalCoordinateNewton local out of +bounds: " << local_iterative;
//                              local_iterative[i] = 0.99;
                        return 1;
 
                }
        }
 
        Verbose(0) << "getLocalCoordinateNewton found local: " << local;
 
        return 0;
}
 
 
template<class T>
std::pair<bool, T> getInterpolatedVertexField( const Eagle::PhysicalSpace::point& world_pos, const Field& vertex_data_field )
{
#ifdef VERBOSE
        std::cout << "getInterpolatedVertexField() ENTER" << std::endl;
#endif
Eagle::PhysicalSpace::point local_pos;
FEMIndex_t cell_nr = 0;
 
//int check = getLocalCoordinate( local_pos, cell_nr, world_pos );
//      Verbose(0) << "getInterpolatedVertexField() local projective: " << local_pos;
 
int check = 1;
 
// get candidates:
std::set<index_t> cell_candidates;
Eagle::PhysicalSpace::tvector range;
        range.x() = 3;
        range.y() = 3;
        range.z() = 3;
 
        getCellCandidates(cell_candidates, range, world_pos );
 
unsigned count = 0;
std::map<double, index_t> closest;
std::map<index_t, Eagle::PhysicalSpace::point> local;
 
        std::set<index_t>::const_iterator cand_it;
        for( cand_it = cell_candidates.begin(); cand_it != cell_candidates.end(); cand_it++ )
        {
                check = getLocalCoordinateNewton( local_pos, *cand_it, world_pos );
                Verbose(0) << "getInterpolatedVertexField() local newtonr: " << local_pos;
 
        double greatest_value = std::max( std::max( local_pos[0], local_pos[1] ), local_pos[2]);
                closest[ greatest_value ] = *cand_it;
                local[*cand_it] = local_pos;
 
                if( check == 0)
                        break;
        }
 
        if( check != 0 )
        {
                Verbose(0)  << "getInterpolatedVertexField()::getLocalCoordinate() returned 1!";
 
                if( closest.size() == 0 )
                        return std::pair<bool, T>(false, T());
 
                cell_nr   = closest.begin()->second;
                local_pos = local[cell_nr];
 
                for( unsigned i = 0; i < 3; i++ )
                {
                        if( local_pos[i] < -1 )
                                local_pos[i] = -0.99;
                        if( local_pos[i] > 1. )
                                local_pos[i] = 0.99;
                }
 
                RefPtr< MemArray<1, T> > data_arr;
 
#ifdef  USE_OPENMP
#pragma omp critical
{
#endif
//      std::cout << "creating averaged field" << std::endl;
//      exit(0);
        data_arr = vertex_data_field.getData();
        assert(data_arr);
#ifdef  USE_OPENMP
}
#endif
 
Verbose(0) << "-------------------------------------------------------------------------------- cell_nr:   " << cell_nr;
Verbose(0) << "-------------------------------------------------------------------------------- local_pos: " << local_pos;
 
                return std::pair<bool, T>(true, FEMFunc::getInterpolatedC3D8<T>(local_pos, *v_cell_arr, cell_nr, *data_arr ));
        }
//      assert(check == 0);
 
#ifdef VERBOSE
        std::cout << "getInterpolatedVertexField()[" << check << "] world: " << world_pos << "-->" << local_pos << std::endl;
#endif
 
//RefPtr<       MemArray<2, FEMIndex_t> > cell_arr = v_cell_positions->getData();
//        assert(v_cell_arr);
//#ifdef        USE_OPENMP
 
//#endif
        RefPtr< MemArray<1, T> > data_arr;
#ifdef  USE_OPENMP
#pragma omp critical
{
#endif
//      std::cout << "creating averaged field" << std::endl;
//      exit(0);
        data_arr = vertex_data_field.getData();
        assert(data_arr);
#ifdef  USE_OPENMP
}
#endif
 
return std::pair<bool, T>(true, FEMFunc::getInterpolatedC3D8<T>( local_pos, *v_cell_arr, cell_nr, *data_arr ));
 
 
}
 
};
 
 
 
//}
 
}
 
 
#endif