Author: cosurgi Date: 2009-03-16 01:20:35 +0100 (Mon, 16 Mar 2009) New Revision: 1720
Added: trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp trunk/pkg/snow/Engine/ElawSnowLayersDeformation.hpp Log: I forgot to add those, sorry. (I wonder how it could stay unnoticed for so long ;) Added: trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp =================================================================== --- trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp 2009-03-13 09:08:31 UTC (rev 1719) +++ trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp 2009-03-16 00:20:35 UTC (rev 1720) @@ -0,0 +1,187 @@ +/************************************************************************* +* Copyright (C) 2009 by Janek Kozicki * +* [email protected] * +* * +* This program is free software; it is licensed under the terms of the * +* GNU General Public License v2 or later. See file LICENSE for details. * +*************************************************************************/ + +#include"ElawSnowLayersDeformation.hpp" +#include<yade/pkg-dem/CohesiveFrictionalBodyParameters.hpp> +#include<yade/pkg-dem/SpheresContactGeometry.hpp> +#include<yade/pkg-dem/CohesiveFrictionalContactInteraction.hpp> +#include<yade/pkg-dem/SDECLinkPhysics.hpp> +#include<yade/core/Omega.hpp> +#include<yade/core/MetaBody.hpp> +#include<yade/pkg-common/Force.hpp> +#include<yade/pkg-common/Momentum.hpp> +#include<yade/core/PhysicalAction.hpp> +#include<yade/pkg-snow/BssSnowGrain.hpp> +#include<yade/pkg-snow/BshSnowGrain.hpp> + + +ElawSnowLayersDeformation::ElawSnowLayersDeformation() : InteractionSolver() , actionForce(new Force) , actionMomentum(new Momentum) +{ + sdecGroupMask=1; + creep_viscosity = 1000.0; +} + + +void ElawSnowLayersDeformation::registerAttributes() +{ + InteractionSolver::registerAttributes(); + REGISTER_ATTRIBUTE(sdecGroupMask); + REGISTER_ATTRIBUTE(creep_viscosity); +} + +void ElawSnowLayersDeformation::action(MetaBody* ncb) +{ + //return; + shared_ptr<BodyContainer>& bodies = ncb->bodies; + +// Real dt = Omega::instance().getTimeStep(); + InteractionContainer::iterator ii = ncb->transientInteractions->begin(); + InteractionContainer::iterator iiEnd = ncb->transientInteractions->end(); + for ( ; ii!=iiEnd ; ++ii ) + { + if ((*ii)->isReal) + { + const shared_ptr<Interaction>& contact = *ii; + int id1 = contact->getId1(); + int id2 = contact->getId2(); + + + +// if(!(((id1 == 17 && id2 == 13) || (id1 == 13 && id2 == 17)))) +// continue; + + if ( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) ) + continue; // skip other groups, + + std::cerr << __FILE__ << " " << id1 << " " << id2 << "\n"; + + CohesiveFrictionalBodyParameters* de1 = YADE_CAST<CohesiveFrictionalBodyParameters*>((*bodies)[id1]->physicalParameters.get()); + CohesiveFrictionalBodyParameters* de2 = YADE_CAST<CohesiveFrictionalBodyParameters*>((*bodies)[id2]->physicalParameters.get()); +// SpheresContactGeometry* currentContactGeometry = YADE_CAST<SpheresContactGeometry*>(contact->interactionGeometry.get()); + CohesiveFrictionalContactInteraction* currentContactPhysics = YADE_CAST<CohesiveFrictionalContactInteraction*> (contact->interactionPhysics.get()); + + BssSnowGrain* b1 = dynamic_cast<BssSnowGrain*>((*bodies)[id1]->interactingGeometry.get()); + BssSnowGrain* b2 = dynamic_cast<BssSnowGrain*>((*bodies)[id2]->interactingGeometry.get()); + + BshSnowGrain* B1 = dynamic_cast<BshSnowGrain*>((*bodies)[id1]->geometricalModel.get()); + BshSnowGrain* B2 = dynamic_cast<BshSnowGrain*>((*bodies)[id2]->geometricalModel.get()); + + Vector3r F = currentContactPhysics->shearForce + currentContactPhysics->normalForce; + + //FIXME: moment is still not used... + Vector3r M = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending; + + if(b1) + { + Vector3r c1 = b1->m_copy.c_axis; + c1 /= c1.Length(); + //FIXME: make sure that F is always in correct direction + F *= -1.0; + Vector3r F1 = de1->se3.orientation.Conjugate() * F; + F1 = F1 - c1*F1.Dot(c1); + b1->m_copy.has_deformed(); + B1->has_deformed(); + + std::vector<Real> layer_totals; + Real all_layers_total(0); + layer_totals.resize(b1->m_copy.slices.size() , 0); + + BOOST_FOREACH(const depth_one& t,b1->depths[id2]) + { + int i = t.i; + //int j = t.j; + + layer_totals[i] += 1.0;//+= t.current_depth - t.original_depth; + ++all_layers_total; + //b1->m_copy.slices[i][j] += F1/(creep_viscosity); + //B1-> slices[i][j] += F1/(creep_viscosity); + } + + for(size_t i = 0 ; i < b1->m_copy.slices.size() ; ++i ) + { + if(layer_totals[i] > 0) + { + for(size_t j = 0 ; j < b1->m_copy.slices[0].size() ; ++j ) + { + b1->m_copy.slices[i][j] += (F1*layer_totals[i])/(creep_viscosity*all_layers_total); + B1-> slices[i][j] += (F1*layer_totals[i])/(creep_viscosity*all_layers_total); + } + } + } + + // for(size_t i=0;i < b1->m_copy.slices.size();++i) + // { + // for(size_t j=0 ; j < b1->m_copy.slices[i].size() ; ++j) + // { + // //if(b1->m_copy.slices[i][j][2]>0) + // { + // b1->m_copy.slices[i][j] += F1/(creep_viscosity); + // B1-> slices[i][j] += F1/(creep_viscosity); + // } + // } + // } + } + + if(b2) + { + Vector3r c2 = b2->m_copy.c_axis; + c2 /= c2.Length(); + F *= -1.0; + Vector3r F2 = de2->se3.orientation.Conjugate() * F; + F2 = F2 - c2*F2.Dot(c2); + b2->m_copy.has_deformed(); + B2->has_deformed(); + + std::vector<Real> layer_totals; + Real all_layers_total(0); + layer_totals.resize(b2->m_copy.slices.size() , 0); + + BOOST_FOREACH(const depth_one& t,b2->depths[id1]) + { + int i = t.i; + //int j = t.j; + + layer_totals[i] += 1.0;//+= t.current_depth - t.original_depth; + // FIXME: divide force between all layers proportionally to sum of (t.current_depth - t.original_depth) of all nodes in each layer + ++all_layers_total; + //b2->m_copy.slices[i][j] += F2/(creep_viscosity); + //B2-> slices[i][j] += F2/(creep_viscosity); + } + + for(size_t i = 0 ; i < b2->m_copy.slices.size() ; ++i ) + { + if(layer_totals[i] > 0) + { + for(size_t j = 0 ; j < b2->m_copy.slices[0].size() ; ++j ) + { + b2->m_copy.slices[i][j] += (F2*layer_totals[i])/(creep_viscosity*all_layers_total); + B2-> slices[i][j] += (F2*layer_totals[i])/(creep_viscosity*all_layers_total); + } + } + } + + + // for(size_t i=0;i < b2->m_copy.slices.size();++i) + // { + // for(size_t j=0 ; j < b2->m_copy.slices[i].size() ; ++j) + // { + // //if(b2->m_copy.slices[i][j][2]>0) + // { + // b2->m_copy.slices[i][j] += F2/(creep_viscosity); + // B2-> slices[i][j] += F2/(creep_viscosity); + // } + // } + // } + } + + } + } +} + +YADE_PLUGIN(); + Added: trunk/pkg/snow/Engine/ElawSnowLayersDeformation.hpp =================================================================== --- trunk/pkg/snow/Engine/ElawSnowLayersDeformation.hpp 2009-03-13 09:08:31 UTC (rev 1719) +++ trunk/pkg/snow/Engine/ElawSnowLayersDeformation.hpp 2009-03-16 00:20:35 UTC (rev 1720) @@ -0,0 +1,41 @@ +/************************************************************************* +* Copyright (C) 2009 by Janek Kozicki * +* [email protected] * +* * +* This program is free software; it is licensed under the terms of the * +* GNU General Public License v2 or later. See file LICENSE for details. * +*************************************************************************/ + +#pragma once + +#include<yade/core/InteractionSolver.hpp> + +#include <set> +#include <boost/tuple/tuple.hpp> + +class PhysicalAction; + +class ElawSnowLayersDeformation : public InteractionSolver +{ +/// Attributes + private : + shared_ptr<PhysicalAction> actionForce; + shared_ptr<PhysicalAction> actionMomentum; + + public : + int sdecGroupMask; + Real creep_viscosity; + + ElawSnowLayersDeformation(); + void action(MetaBody*); + + protected : + void registerAttributes(); + NEEDS_BEX("Force","Momentum"); + REGISTER_CLASS_NAME(ElawSnowLayersDeformation); + REGISTER_BASE_CLASS_NAME(InteractionSolver); +}; + +REGISTER_SERIALIZABLE(ElawSnowLayersDeformation); + + _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
