FEMFields.hpp

 
//
// $Id: hpp,v 1.0 2000/01/01 00:00:01 werner Exp $
//
// $Log: hpp,v $
// Created 2012/05/21 12:57:24  marcel
// Initial version
//
//
#ifndef __FEMFIELDS_HPP
#define __FEMFIELDS_HPP
 
#include "FEM.hpp"
#include <eagle/PhysicalSpace.hpp>
#include "gridtypesDllApi.h"
#include <grid/Grid.hpp>
 
//#include <ocean/plankton/VInput.hpp>
 
namespace Fiber
{
 
struct FEMOperatorData
{
        WeakPtr<Grid>           m_vgrid;
        WeakPtr<Grid>           m_igrid;
 
        string m_field_name;
        size_t m_size; 
 
        FEMOperatorData( WeakPtr<Grid> vgrid, WeakPtr<Grid> igrid, string field_name, size_t size = 0 )
        : m_vgrid( vgrid )
        , m_igrid( igrid )
        , m_field_name( field_name )
        , m_size( size )
        {}
};
 
 
struct gridtypes_API IntegrationPointVertices
{
typedef MemArray<1, Eagle::PhysicalSpace::point> ResultArray_t; 
typedef FEMOperatorData                          Constructor_t;
 
        RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > m_data; 
 
        IntegrationPointVertices(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~IntegrationPointVertices() {} 
 
        RefPtr<ResultArray_t> result() const;
};
 
// stress related fields:       
struct gridtypes_API MeanStress
{
typedef MemArray<1, float>      ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<MemArray<1, float> > m_data;
 
        MeanStress(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~MeanStress() {}
 
        RefPtr<ResultArray_t> result() const;
};
 
 
struct gridtypes_API Invariant1
{
typedef MemArray<1, float>      ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<MemArray<1, float> > m_data;
 
        Invariant1(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~Invariant1() {}
 
        RefPtr<ResultArray_t> result() const;
};
 
struct gridtypes_API Invariant2
{
typedef MemArray<1, float>      ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<MemArray<1, float> > m_data;
 
        Invariant2(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~Invariant2() {}
 
        RefPtr<ResultArray_t> result() const;
};
 
struct gridtypes_API Invariant3
{
typedef MemArray<1, float>      ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<MemArray<1, float> > m_data;
 
        Invariant3(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~Invariant3() {}
 
        RefPtr<ResultArray_t> result() const;
};
 
struct gridtypes_API Mises
{
typedef MemArray<1, float>      ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<MemArray<1, float> > m_data;
 
        Mises(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~Mises() {}
 
        RefPtr<ResultArray_t> result() const;
};
 
struct gridtypes_API DeviatoricStress
{
typedef MemArray<1, Eagle::metric33>    ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<MemArray<1, Eagle::metric33> > m_data;
 
        DeviatoricStress(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~DeviatoricStress() {}
 
        RefPtr<ResultArray_t> result() const;
};
 
struct gridtypes_API DeviatoricRadius
{
typedef MemArray<1, float>      ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<MemArray<1, float> > m_data;
 
        DeviatoricRadius(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
        ~DeviatoricRadius() {}
 
        RefPtr<ResultArray_t> result() const;
};
 
template<class T>  void initValues( T& foo);
template<> void initValues<float>( float& foo);
template<> void initValues<double>( double& foo);
template<> void initValues<Eagle::PhysicalSpace::tvector>( Eagle::PhysicalSpace::tvector& foo);
template<> void initValues<Eagle::metric33>( Eagle::metric33& foo);
 
 
template<class T>
struct gridtypes_API AverageField
{
typedef MemArray<1, T>     ResultArray_t; 
typedef FEMOperatorData    Constructor_t; 
 
RefPtr<ResultArray_t> m_data;
 
        AverageField(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() )
        {
                std::cout << "AverageStress::AverageStress()" << std::endl; 
                std::cout.flush(); 
 
                RefPtr<Skeleton> v_cells =    (*FEMData.m_vgrid)( FEM::cellID3D() ); 
                RefPtr<Skeleton> v_vertices = FEMData.m_vgrid->findVertices(); 
                assert(v_cells); 
                assert(v_vertices); 
 
                RefPtr<Representation> v_cells_as_vertices = (*v_cells)(v_vertices); 
                assert(v_cells_as_vertices); 
 
//RefPtr<Representation> v_vertices_as_cells = (*v_vertices)(v_cells); 
//        assert(v_vertices_as_cells); 
 
                RefPtr<Field> v_c_v_pos = v_cells_as_vertices->Positions(); 
                assert(v_c_v_pos); 
 
//RefPtr<Field> v_v_c_pos = v_vertices_as_cells->Positions(); 
//        assert(v_v_c_pos); 
 
                RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > vert_arr = FEMData.m_vgrid->getCartesianPositions()->getData(); 
                assert(vert_arr); 
 
                RefPtr<Field> i_field = (*FEMData.m_igrid)( FEMData.m_field_name ); 
                assert(i_field); 
 
                RefPtr<MemArray<1, T> > i_av_source_arr = i_field->getData(); 
                assert(i_av_source_arr); 
                MultiArray<1, T>&i_av_values = *i_av_source_arr; 
 
                m_data = new MemArray<1, T>(vert_arr->Size() ); 
                MultiArray<1, T>& av_values = *m_data; 
                MultiIndex<1> mi; 
                mi = 0; 
                do
                {
                        //Wizt::VInputValueTrait<T>::init( av_values[mi] );
                        initValues<T>( av_values[mi] );
 
                }while( mi.inc( vert_arr->Size() ) ); 
 
                RefPtr<Skeleton> i_vertices =  FEMData.m_igrid->findVertices(); 
                RefPtr<Skeleton> i_cells = (*FEMData.m_igrid)( FEM::cellID3D() ); 
                assert(i_vertices); 
                assert(i_cells); 
                RefPtr<Representation> i_cells_as_vertices = (*i_cells)(i_vertices); 
                assert(i_cells_as_vertices); 
                RefPtr<Field> i_c_v_pos = i_cells_as_vertices->Positions(); 
                assert(i_c_v_pos); 
 
                {
                        RefPtr<MemArray<1, std::vector<uint32_t> > > cell_v_arr = v_c_v_pos->getData(); 
                        RefPtr<MemArray<1, std::vector<uint32_t> > > cell_i_arr = i_c_v_pos->getData(); 
 
                        if( cell_v_arr && cell_i_arr )
                        {
                                MultiArray<1, std::vector<uint32_t> >& cell_v = *cell_v_arr; 
                                MultiArray<1, std::vector<uint32_t> >& cell_i = *cell_i_arr; 
 
// iterate over cells 
                                mi = 0; 
                                do
                                {
                                        std::vector<uint32_t>& one_vcell = cell_v[mi]; 
                                        std::vector<uint32_t>& one_icell = cell_i[mi]; 
 
                                        // sum up: 
                                        if( one_icell.size() == 1 )  // reduced integrated cells
                                        {
                                                // for vertices of cells compute contribution of reduced stress value 
                                                // to each vertex equally
                                                for( size_t i = 0; i < one_vcell.size(); i++ )
                                                {
                                                const unsigned& mv = one_vcell[i]; 
                                                const unsigned& mj = one_icell[0]; 
 
                                                        av_values[mv] += i_av_values[mj] / one_vcell.size();
                                                }
                                        } 
//              else 
                                        if( one_icell.size() == 27 && one_vcell.size() == 20 ) // c3d20 element 
                                        {
                                                for( size_t i = 0; i < one_vcell.size(); i++ )
                                                {
                                                const unsigned& mv = one_vcell[i]; 
                                                double weight_sum = 0; 
 
                                                T tmp; 
                                                        initValues<T>(tmp);
 
                                                        for( size_t j = 0; j < one_icell.size(); j++ )
                                                        {
                                                        const unsigned& mj = one_icell[j]; 
 
                                                        double weight = FEMFunc::getC3D20GaussWeight(j, i); 
 
                                                                tmp += T(i_av_values[mj] * weight); 
 
                                                                weight_sum += weight; 
                                                        } 
                                                        av_values[mv] +=  tmp / weight_sum;
                                                }
                                        } 
                                        if( one_icell.size() == 8 && one_vcell.size() == 8 ) // c3d8 element
                                        {
                                                for( size_t i = 0; i < one_vcell.size(); i++ )
                                                {
                                                const unsigned& mv = one_vcell[i]; 
 
                                                double weight_sum = 0; 
 
                                                T tmp; 
                                                        initValues<T>(tmp);
 
                                                        for( size_t j = 0; j < one_icell.size(); j++ )
                                                        {
                                                        const unsigned& mj = one_icell[j]; 
 
                                                        double weight = FEMFunc::getC3D8GaussWeight(j, i); 
                                                                tmp += T(i_av_values[mj] * weight); 
                                                                weight_sum += weight; 
                                                        } 
                                                        av_values[mv] += tmp/weight_sum; 
                                
                                                }
                                        } 
 
                                }while( mi.inc( cell_v_arr->Size() ) ); 
 
                                std::cout << "AverageStress::AverageStress()END" << std::endl; 
                                std::cout.flush();
                        } 
                } 
                // fixed sizes? 
                {
                        RefPtr<MemArray<2, uint32_t> > cell_v_arr = v_c_v_pos->getData();; 
                        RefPtr<MemArray<2, uint32_t> > cell_i_arr = i_c_v_pos->getData();; 
                        if( cell_v_arr && cell_i_arr )
                        {
                                MultiArray<2, uint32_t>& coords_cell     = *cell_v_arr; 
                                MultiArray<2, uint32_t>& int_points_cell = *cell_i_arr; 
 
                                MultiIndex<2> mi, mv; 
                                mi = 0, 0; 
                                mv = 0, 0; 
 
                                std::vector<double> weight_sum( vert_arr->Size()[0], 0.0 ); 
 
                                // weighted sum: 
                                // for cells 
                                for( index_t c = 0; c < cell_v_arr->Size()[1]; c++)
                                {
                                        // for vertices (of cell)
                                        for( index_t v = 0; v < cell_v_arr->Size()[0]; v++ )
                                        {
                                                mv = v, c; 
 
                                                // for int Eagle::points 
                                                for( index_t i = 0; i < cell_i_arr->Size()[0]; i++ )
                                                {
                                                        mi = i, c; 
                                                double weight = 0; 
 
                                                        if( cell_v_arr->Size()[0] == 20 && cell_i_arr->Size()[0] == 27 )
                                                                weight = FEMFunc::getC3D20GaussWeight(i, v); 
                                                        else if( cell_v_arr->Size()[0] == 8 && cell_i_arr->Size()[0] == 8 )
                                                                weight = FEMFunc::getC3D8GaussWeight(i, v); 
                                                        else if( cell_v_arr->Size()[0] == 8 && cell_i_arr->Size()[0] == 1  )
                                                                weight = 1; 
                                                        else if( cell_v_arr->Size()[0] == 20 && cell_i_arr->Size()[0] == 8  )
                                                                weight = FEMFunc::getC3D20RGaussWeight(i, v); 
                                                        else if( cell_v_arr->Size()[0] == 4 && cell_i_arr->Size()[0] == 1  )
                                                                weight = 1; 
                                                        else
                                                                assert( 0 && "unsupported element type" ); 
 
                                                        av_values[ coords_cell[mv] ]  += T(i_av_values[ int_points_cell[mi] ] * weight); 
                                                        weight_sum[ coords_cell[mv] ] += weight; 
                                                }
                                        }
                                } 
 
                                // normalize: 
                                for( index_t i = 0; i < vert_arr->Size()[0]; i++ )
                                        av_values[ i ] /= weight_sum[ i ]; 
 
                        }
                }
        }
        
 
        ~AverageField() {}
 
        RefPtr<ResultArray_t> result() const
        {
                return m_data;
        }
};
 
 
template<class T>
struct gridtypes_API ComputeCellField
{
typedef MemArray<2, T>     ResultArray_t; 
typedef FEMOperatorData    Constructor_t; 
 
RefPtr<ResultArray_t> m_data;
 
        ComputeCellField(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() )
        {
#ifdef VERBOSE
                std::cout << "ComputeCellField::ComputeCellField()" << std::endl; 
                std::cout.flush(); 
#endif
//              std::cin.ignore();
                
                RefPtr<Skeleton> v_cells =    (*FEMData.m_vgrid)( FEM::cellID3D() ); 
                RefPtr<Skeleton> v_vertices = FEMData.m_vgrid->findVertices(); 
                assert(v_cells); 
                assert(v_vertices); 
 
                RefPtr<Representation> v_cells_as_vertices = (*v_cells)(v_vertices); 
                assert(v_cells_as_vertices); 
 
//RefPtr<Representation> v_vertices_as_cells = (*v_vertices)(v_cells); 
//        assert(v_vertices_as_cells); 
 
                RefPtr<Field> v_c_v_pos = v_cells_as_vertices->Positions(); 
                assert(v_c_v_pos); 
 
//RefPtr<Field> v_v_c_pos = v_vertices_as_cells->Positions(); 
//        assert(v_v_c_pos); 
 
                RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > vert_arr = FEMData.m_vgrid->getCartesianPositions()->getData(); 
                assert(vert_arr); 
 
                RefPtr<Field> i_field = (*FEMData.m_igrid)( FEMData.m_field_name ); 
                assert(i_field); 
 
                RefPtr<MemArray<1, T> > i_av_source_arr = i_field->getData(); 
                assert(i_av_source_arr); 
                MultiArray<1, T>&i_av_values = *i_av_source_arr; 
 
                RefPtr<MemArray<2, uint32_t> > cell_v_arr = v_c_v_pos->getData(); 
                assert( cell_v_arr );
 
                m_data = new MemArray<2, T>( cell_v_arr->Size() ); 
                MultiArray<2, T>& new_values = *m_data; 
                MultiIndex<2> mi; 
                mi = 0,0; 
                do
                {
                        //Wizt::VInputValueTrait<T>::init( av_values[mi] );
                        initValues<T>( new_values[mi] );
 
                }while( mi.inc( cell_v_arr->Size() ) ); 
 
                RefPtr<Skeleton> i_vertices =  FEMData.m_igrid->findVertices(); 
                RefPtr<Skeleton> i_cells = (*FEMData.m_igrid)( FEM::cellID3D() ); 
                assert(i_vertices); 
                assert(i_cells); 
                RefPtr<Representation> i_cells_as_vertices = (*i_cells)(i_vertices); 
                assert(i_cells_as_vertices); 
                RefPtr<Field> i_c_v_pos = i_cells_as_vertices->Positions(); 
                assert(i_c_v_pos); 
 
// fixed size support only (nor now)
                {
//                      RefPtr<MemArray<2, uint32_t> > cell_v_arr = v_c_v_pos->getData(); 
                        RefPtr<MemArray<2, uint32_t> > cell_i_arr = i_c_v_pos->getData(); 
                        if( cell_v_arr && cell_i_arr )
                        {
//                              MultiArray<2, uint32_t>& coords_cell     = *cell_v_arr; 
                                MultiArray<2, uint32_t>& int_points_cell = *cell_i_arr; 
 
                                MultiIndex<2> mi, mv; 
                                mi = 0, 0; 
                                mv = 0, 0; 
 
//                              vector<double> weight_sum( vert_arr->Size()[0], 0.0 ); 
                                RefPtr<MemArray<2, double> > weight_sum_arr = new MemArray<2, double> ( cell_v_arr->Size() ); 
                                MultiArray<2, double>& weight_sum = *weight_sum_arr;
 
                                // weighted sum: 
                                // for cells 
                                for( index_t c = 0; c < cell_v_arr->Size()[1]; c++)
                                {
                                        // for vertices (of cell)
                                        for( index_t v = 0; v < cell_v_arr->Size()[0]; v++ )
                                        {
                                                mv = v, c; 
 
                                                // for int Eagle::points 
                                                for( index_t i = 0; i < cell_i_arr->Size()[0]; i++ )
                                                {
                                                        mi = i, c; 
                                                double weight = 0; 
 
                                                        if( cell_v_arr->Size()[0] == 20 && cell_i_arr->Size()[0] == 27 )
                                                                weight = FEMFunc::getC3D20GaussWeight(i, v); 
                                                        else if( cell_v_arr->Size()[0] == 8 && cell_i_arr->Size()[0] == 8 )
                                                                weight = FEMFunc::getC3D8GaussWeight(i, v); 
                                                        else if( cell_v_arr->Size()[0] == 8 && cell_i_arr->Size()[0] == 1  )
                                                                weight = 1; 
                                                        else if( cell_v_arr->Size()[0] == 20 && cell_i_arr->Size()[0] == 8  )
                                                                weight = FEMFunc::getC3D20RGaussWeight(i, v); 
                                                        else if( cell_v_arr->Size()[0] == 4 && cell_i_arr->Size()[0] == 1  )
                                                                weight = 1; 
                                                        else
                                                                assert( 0 && "unsupported element type" ); 
 
                                                        new_values[ mv ] += T(i_av_values[ int_points_cell[mi] ] * weight); 
                                                        weight_sum[ mv ] += weight; 
                                                }
                                        }
                                } 
 
                                std::cout << "ComputeCellField:: normalize" << std::endl; 
                                // normalize
                                mv = 0,0; 
                                do
                                {
                                        std::cout << mv << std::endl;
                                        new_values[ mv ] /= weight_sum[mv];
                                }while( mv.inc( cell_v_arr->Size() ) ); 
 
                                     }
                        }
                }
        
 
        ~ComputeCellField() {}
 
        RefPtr<ResultArray_t> result() const
        {
                return m_data;
        }
};
 
struct gridtypes_API VertexCellNeighborhood
{
typedef MemArray<1, FEMIndices_t>       ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<ResultArray_t> m_data;
 
        VertexCellNeighborhood(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
//      ~CellNeighborhood() {} 
 
RefPtr<ResultArray_t> result() const;
};
 
 
struct gridtypes_API CellNeighborhoodVariable
{
typedef MemArray<1, FEMIndices_t>       ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<ResultArray_t> m_data;
 
        CellNeighborhoodVariable(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
//      ~CellNeighborhood() {} 
 
RefPtr<ResultArray_t> result() const;
};
 
struct gridtypes_API CellDirectionsVariable
{
typedef MemArray<1, FEMDirections_t>    ResultArray_t; 
typedef FEMOperatorData         Constructor_t; 
 
RefPtr<ResultArray_t> m_data;
 
        CellDirectionsVariable(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
//      ~CellNeighborhood() {} 
 
RefPtr<ResultArray_t> result() const;
};
 
 
struct gridtypes_API CellNeighborhoodFixed
{
typedef MemArray<2, FEMIndex_t>  ResultArray_t; 
typedef FEMOperatorData          Constructor_t; 
 
RefPtr<ResultArray_t> m_data;
 
        CellNeighborhoodFixed(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
//      ~CellNeighborhood() {} 
 
RefPtr<ResultArray_t> result() const;
};
 
struct gridtypes_API CellDirectionsFixed
{
typedef MemArray<2, Eagle::PhysicalSpace::tvector> ResultArray_t; 
typedef FEMOperatorData                            Constructor_t; 
 
RefPtr<ResultArray_t> m_data;
 
        CellDirectionsFixed(const FEMOperatorData& FEMData, const MemBase::Creator_t&C = MemCore::NullPtr() ); 
//      ~CellNeighborhood() {} 
 
RefPtr<ResultArray_t> result() const;
};
 
}
 
#endif