yade-dev team mailing list archive

[chareyre at BerliOS <chareyre@xxxxxxxxxxxxxxx>]

Author: chareyre
Date: 2009-07-24 00:58:43 +0200 (Fri, 24 Jul 2009)
New Revision: 1881

Modified:
   trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.cpp
   trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.hpp
Log:
Update of the read/write functions for the connection 
with comsol. (coupling via Berkeley<->Grenoble mail servers! ;)) 



Modified: trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.cpp	2009-07-21 16:04:11 UTC (rev 1880)
+++ trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.cpp	2009-07-23 22:58:43 UTC (rev 1881)
@@ -1,6 +1,6 @@
 /*************************************************************************
-*  Copyright (C) 2004 by Janek Kozicki                                   *
-*  cosurgi@xxxxxxxxxx                                                    *
+*  Copyright (C) 2004 by Andrea Cortis & Bruno Chareyre                  *
+*  acortis@xxxxxxx,   bruno.chareyre@xxxxxxxxxxx                         *
 *                                                                        *
 *  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. *
@@ -8,17 +8,24 @@
 
 #include "HydraulicForceEngine.hpp"
 #include <yade/pkg-common/ParticleParameters.hpp>
-#include<yade/core/MetaBody.hpp>
+#include <yade/core/MetaBody.hpp>
+#include <yade/core/Body.hpp>
 #include <yade/pkg-dem/CohesiveFrictionalBodyParameters.hpp>
 #include <vector>
+#include "HydraulicForceEngine.hpp"
+#include <yade/pkg-common/ParticleParameters.hpp>
+#include <yade/core/MetaBody.hpp>
+#include <yade/pkg-dem/CohesiveFrictionalBodyParameters.hpp>
+#include <vector>
+#include <yade/extra/Shop.hpp>
+#include<yade/pkg-common/InteractingSphere.hpp>
 
-long int HFinversion_counter = 0;
-bool HFinverted = false;
 vector<Real> initialPositions;
 
 HydraulicForceEngine::HydraulicForceEngine() : gravity(Vector3r::ZERO), isActivated(false)
 {
 dummyParameter = false;
+
 }
 
 
@@ -35,22 +42,24 @@
     {
         shared_ptr<BodyContainer>& bodies = ncb->bodies;
         
+        
+        ///This commented is to artificially remove eroded (i.e. bond-breaked) grains by applying a force moving them away
 //cerr << "HFinverted " << HFinverted << " Omega::instance().getTimeStep() " << Omega::instance().getCurrentIteration() << endl;
-        if (HFinversion_counter < (Omega::instance().getCurrentIteration() - 1500))
-        {
-        //cerr << "HFinverted " << HFinverted << " HFinversion_counter " << HFinversion_counter << endl;
-            HFinversion_counter = Omega::instance().getCurrentIteration();
-            HFinverted = !HFinverted;
-//             initialPositions.clear();
-//             for (int i =0; i < initialPositions.size(); ++i)
-//             {
-//             	initialPositions[i-6] = (static_cast<CohesiveFrictionalBodyParameters*> (b->physicalParameters.get()))->se3.position[1];
-//             	
-//             }
-            
-            
-        }
-        
+//         if (HFinversion_counter < (Omega::instance().getCurrentIteration() - 1500))
+//         {
+//         //cerr << "HFinverted " << HFinverted << " HFinversion_counter " << HFinversion_counter << endl;
+//             HFinversion_counter = Omega::instance().getCurrentIteration();
+//             HFinverted = !HFinverted;
+// //             initialPositions.clear();
+// //             for (int i =0; i < initialPositions.size(); ++i)
+// //             {
+// //             	initialPositions[i-6] = (static_cast<CohesiveFrictionalBodyParameters*> (b->physicalParameters.get()))->se3.position[1];
+// //             	
+// //             }
+//             
+//             
+//         }
+//         
 
         if (HFinverted)
         {
@@ -74,29 +83,62 @@
         {
             long Nspheres;
             long id;
+	    Real temp;
+	    Real realID;
             Real fx, fy, fz, mx, my, mz;
-            ifstream file ("/home/bruno/YADE/data/hydraulic_actions.dat");
-            file >> Nspheres;
-            initialPositions.resize(Nspheres);
+            ifstream file (inputFile.c_str());
+	    ofstream ofile (outputFile.c_str());
+            
             //cerr << "r = " << R << endl;
+		#define USELINENUMBERSFORID
             if (file.is_open())
             {
-                for (long i=1; i<Nspheres; ++i)
-                {
+		file >> realID >> temp >> temp>> temp>> temp>> temp>> temp;
+		Nspheres = realID; 
+		cerr << "Nspheres" << Nspheres << endl;
+		initialPositions.resize(Nspheres);
+		
+		for (long i=6; i<Nspheres+6; ++i)
+		{
+			file >> realID >> fx >> fy >> fz >> mx >> my >> mz;
+			id = realID;
+		#ifdef USELINENUMBERSFORID
+                	id = i;
+		#endif
+                
 
-                    file >> id >> fx >> fy >> fz >> mx >> my >> mz;
-
                     Vector3r f (fx,fy,fz);
                     Vector3r t (mx,my,mz);
+		    f*=forceFactor;
+		    t*=forceFactor;
+// 		    cerr << "f=" << f << " on " << id << endl;
                     //f /= -10000;
                     //t *= 0;
 						  	ncb->bex.addForce(id,f);
 							ncb->bex.addTorque(id,t);
-
                 }
+		file.close();
             }
             else cerr << "problem opening the file with hydraulic forces" << endl;
-            file.close();
+            
+	    if (savePositions)
+	    #ifdef USELINENUMBERSFORID
+	    	{Shop::saveSpheresToFile(outputFile.c_str()); savePositions = false;}//only once
+	    #else
+	    {
+		    const shared_ptr<MetaBody>& rootBody=Omega::instance().getRootBody();
+		    ofstream f(outputFile.c_str());
+		    if(!f.good()) throw runtime_error("Unable to open file `"+outputFile+"'");
+		    FOREACH(shared_ptr<Body> b, *rootBody->bodies){
+			    if (!b->isDynamic) continue;
+			    shared_ptr<InteractingSphere>	intSph=dynamic_pointer_cast<InteractingSphere>(b->interactingGeometry);
+			    if(!intSph) continue;
+			    const Vector3r& pos=b->physicalParameters->se3.position;
+			    f<< b->getId()<<" "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<" "<<intSph->radius<<endl; // <<" "<<1<<" "<<1<<endl;
+		    }
+		    f.close();
+	    }
+	    #endif
 
         }
 

Modified: trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.hpp	2009-07-21 16:04:11 UTC (rev 1880)
+++ trunk/pkg/dem/Engine/DeusExMachina/HydraulicForceEngine.hpp	2009-07-23 22:58:43 UTC (rev 1881)
@@ -1,6 +1,6 @@
 /*************************************************************************
-*  Copyright (C) 2004 by Janek Kozicki                                   *
-*  cosurgi@xxxxxxxxxx                                                    *
+*  Copyright (C) 2004 by Andrea Cortis & Bruno Chareyre                  *
+*  acortis@xxxxxxx,   bruno.chareyre@xxxxxxxxxxx                         *
 *                                                                        *
 *  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. *
@@ -16,12 +16,18 @@
 		Vector3r gravity;
 		bool isActivated;
 		bool dummyParameter;
+		bool HFinverted;
+		bool savePositions;
+		string outputFile;
+		string inputFile;
+		int HFinversion_counter;
+		Real forceFactor;
+	
 		HydraulicForceEngine();
 		virtual ~HydraulicForceEngine();
-	
 		virtual void applyCondition(MetaBody*);
 	
-	REGISTER_ATTRIBUTES(DeusExMachina,(gravity)(isActivated)(dummyParameter));
+	REGISTER_ATTRIBUTES(DeusExMachina,(gravity)(isActivated)(dummyParameter)(HFinverted)(savePositions)(outputFile)(inputFile)(HFinversion_counter)(forceFactor));
 	REGISTER_CLASS_NAME(HydraulicForceEngine);
 	REGISTER_BASE_CLASS_NAME(DeusExMachina);
 };