yade-dev team mailing list archive

[cosurgi at BerliOS <cosurgi@xxxxxxxxxxxxxxx>]

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                                   *
+*  cosurgi@xxxxxxxxxx                                                    *
+*                                                                        *
+*  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                                   *
+*  cosurgi@xxxxxxxxxx                                                    *
+*                                                                        *
+*  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);
+
+