doxygen: add newline after \file
[blender.git] / intern / elbeem / intern / ntl_geometryobject.cpp
1 /** \file
2  * \ingroup elbeem
3  */
4 /******************************************************************************
5  *
6  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
7  * Copyright 2003-2006 Nils Thuerey
8  *
9  * a geometry object
10  * all other geometry objects are derived from this one
11  *
12  *****************************************************************************/
13
14
15 #include "ntl_geometryobject.h"
16 #include "ntl_world.h"
17 #include "ntl_matrices.h"
18
19 // for FGI
20 #include "elbeem.h"
21
22 #define TRI_UVOFFSET (1./4.)
23 //#define TRI_UVOFFSET (1./3.)
24
25
26 /*****************************************************************************/
27 /* Default constructor */
28 /*****************************************************************************/
29 ntlGeometryObject::ntlGeometryObject() :
30         mIsInitialized(false), mpMaterial( NULL ),
31         mMaterialName( "default" ),
32         mCastShadows( 1 ), mReceiveShadows( 1 ),
33         mGeoInitType( 0 ), 
34         mInitialVelocity(0.0), mcInitialVelocity(0.0), mLocalCoordInivel(false),
35         mGeoInitIntersect(false),
36         mGeoPartSlipValue(0.0),
37         mcGeoImpactFactor(1.),
38         mVolumeInit(VOLUMEINIT_VOLUME),
39         mInitialPos(0.),
40         mcTrans(0.), mcRot(0.), mcScale(1.),
41         mIsAnimated(false),
42         mMovPoints(), mMovNormals(),
43         mHaveCachedMov(false),
44         mCachedMovPoints(), mCachedMovNormals(),
45         mTriangleDivs1(), mTriangleDivs2(),
46         mMovPntsInited(-100.0), mMaxMovPnt(-1),
47         mcGeoActive(1.),
48         mCpsTimeStart(0.), mCpsTimeEnd(1.0), mCpsQuality(10.),
49         mcAttrFStr(0.),mcAttrFRad(0.), mcVelFStr(0.), mcVelFRad(0.)
50
51 };
52
53
54 /*****************************************************************************/
55 /* Default destructor */
56 /*****************************************************************************/
57 ntlGeometryObject::~ntlGeometryObject() 
58 {
59 }
60
61 /*! is the mesh animated? */
62 bool ntlGeometryObject::getMeshAnimated() {
63         // off by default, on for e.g. ntlGeometryObjModel
64         return false; 
65 }
66
67 /*! init object anim flag */
68 bool ntlGeometryObject::checkIsAnimated() {
69         if(    (mcTrans.accessValues().size()>1)  // VALIDATE
70             || (mcRot.accessValues().size()>1) 
71             || (mcScale.accessValues().size()>1) 
72             || (mcGeoActive.accessValues().size()>1) 
73                         // mcGeoImpactFactor only needed when moving
74             || (mcInitialVelocity.accessValues().size()>1) 
75                 ) {
76                 mIsAnimated = true;
77         }
78
79         // fluid objects always have static init!
80         if(mGeoInitType==FGI_FLUID) {
81                 mIsAnimated=false;
82         }
83         //errMsg("ntlGeometryObject::checkIsAnimated","obj="<<getName()<<" debug: trans:"<<mcTrans.accessValues().size()<<" rot:"<<mcRot.accessValues().size()<<" scale:"<<mcScale.accessValues().size()<<" active:"<<mcGeoActive.accessValues().size()<<" inivel:"<<mcInitialVelocity.accessValues().size()<<". isani?"<<mIsAnimated ); // DEBUG
84         return mIsAnimated;
85 }
86
87 /*****************************************************************************/
88 /* Init attributes etc. of this object */
89 /*****************************************************************************/
90 #define GEOINIT_STRINGS  10
91 static const char *initStringStrs[GEOINIT_STRINGS] = {
92         "fluid",
93         "bnd_no","bnd_noslip",
94         "bnd_free","bnd_freeslip",
95         "bnd_part","bnd_partslip",
96         "inflow", "outflow", "control",
97 };
98 static int initStringTypes[GEOINIT_STRINGS] = {
99         FGI_FLUID,
100         FGI_BNDNO, FGI_BNDNO,
101         FGI_BNDFREE, FGI_BNDFREE,
102         FGI_BNDPART, FGI_BNDPART,
103         FGI_MBNDINFLOW, FGI_MBNDOUTFLOW, 
104         FGI_CONTROL
105 };
106 void ntlGeometryObject::initialize(ntlRenderGlobals *glob) 
107 {
108         //debugOut("ntlGeometryObject::initialize: '"<<getName()<<"' ", 10);
109         // initialize only once...
110         if(mIsInitialized) return;
111         
112         // init material, always necessary
113         searchMaterial( glob->getMaterials() );
114         
115         this->mGeoInitId = mpAttrs->readInt("geoinitid", this->mGeoInitId,"ntlGeometryObject", "mGeoInitId", false);
116         mGeoInitIntersect = mpAttrs->readInt("geoinit_intersect", mGeoInitIntersect,"ntlGeometryObject", "mGeoInitIntersect", false);
117         string ginitStr = mpAttrs->readString("geoinittype", "", "ntlGeometryObject", "mGeoInitType", false);
118         if(this->mGeoInitId>=0) {
119                 bool gotit = false;
120                 for(int i=0; i<GEOINIT_STRINGS; i++) {
121                         if(ginitStr== initStringStrs[i]) {
122                                 gotit = true;
123                                 mGeoInitType = initStringTypes[i];
124                         }
125                 }
126
127                 if(!gotit) {
128                         errFatal("ntlGeometryObject::initialize","Obj '"<<mName<<"', Unknown 'geoinittype' value: '"<< ginitStr <<"' ", SIMWORLD_INITERROR);
129                         return;
130                 }
131         }
132
133         int geoActive = mpAttrs->readInt("geoinitactive", 1,"ntlGeometryObject", "geoActive", false);
134         if(!geoActive) {
135                 // disable geo init again...
136                 this->mGeoInitId = -1;
137         }
138         mInitialVelocity  = vec2G( mpAttrs->readVec3d("initial_velocity", vec2D(mInitialVelocity),"ntlGeometryObject", "mInitialVelocity", false));
139         if(getAttributeList()->exists("initial_velocity") || (!mcInitialVelocity.isInited()) ) {
140                 mcInitialVelocity = mpAttrs->readChannelVec3f("initial_velocity");
141         }
142         // always use channel
143         if(!mcInitialVelocity.isInited()) { mcInitialVelocity = AnimChannel<ntlVec3Gfx>(mInitialVelocity); }
144         mLocalCoordInivel = mpAttrs->readBool("geoinit_localinivel", mLocalCoordInivel,"ntlGeometryObject", "mLocalCoordInivel", false);
145
146         mGeoPartSlipValue = mpAttrs->readFloat("geoinit_partslip", mGeoPartSlipValue,"ntlGeometryObject", "mGeoPartSlipValue", false);
147         bool mOnlyThinInit = false; // deprecated!
148         mOnlyThinInit     = mpAttrs->readBool("geoinit_onlythin", mOnlyThinInit,"ntlGeometryObject", "mOnlyThinInit", false);
149         if(mOnlyThinInit) mVolumeInit = VOLUMEINIT_SHELL;
150         mVolumeInit     = mpAttrs->readInt("geoinit_volumeinit", mVolumeInit,"ntlGeometryObject", "mVolumeInit", false);
151         if((mVolumeInit<VOLUMEINIT_VOLUME)||(mVolumeInit>VOLUMEINIT_BOTH)) mVolumeInit = VOLUMEINIT_VOLUME;
152
153         // moving obs correction factor
154         float impactfactor=1.;
155         impactfactor = (float)mpAttrs->readFloat("impactfactor", impactfactor,"ntlGeometryObject", "impactfactor", false);
156         if(getAttributeList()->exists("impactfactor") || (!mcGeoImpactFactor.isInited()) ) {
157                 mcGeoImpactFactor = mpAttrs->readChannelSinglePrecFloat("impactfactor");
158         }
159
160         // override cfg types
161         mVisible = mpAttrs->readBool("visible", mVisible,"ntlGeometryObject", "mVisible", false);
162         mReceiveShadows = mpAttrs->readBool("recv_shad", mReceiveShadows,"ntlGeometryObject", "mReceiveShadows", false);
163         mCastShadows = mpAttrs->readBool("cast_shad", mCastShadows,"ntlGeometryObject", "mCastShadows", false);
164
165         // read mesh animation channels
166         ntlVec3d translation(0.0);
167         translation = mpAttrs->readVec3d("translation", translation,"ntlGeometryObject", "translation", false);
168         if(getAttributeList()->exists("translation") || (!mcTrans.isInited()) ) {
169                 mcTrans = mpAttrs->readChannelVec3f("translation");
170         }
171         ntlVec3d rotation(0.0);
172         rotation = mpAttrs->readVec3d("rotation", rotation,"ntlGeometryObject", "rotation", false);
173         if(getAttributeList()->exists("rotation") || (!mcRot.isInited()) ) {
174                 mcRot = mpAttrs->readChannelVec3f("rotation");
175         }
176         ntlVec3d scale(1.0);
177         scale = mpAttrs->readVec3d("scale", scale,"ntlGeometryObject", "scale", false);
178         if(getAttributeList()->exists("scale") || (!mcScale.isInited()) ) {
179                 mcScale = mpAttrs->readChannelVec3f("scale");
180         }
181
182         float geoactive=1.;
183         geoactive = (float)mpAttrs->readFloat("geoactive", geoactive,"ntlGeometryObject", "geoactive", false);
184         if(getAttributeList()->exists("geoactive") || (!mcGeoActive.isInited()) ) {
185                 mcGeoActive = mpAttrs->readChannelSinglePrecFloat("geoactive");
186         }
187         // always use channel
188         if(!mcGeoActive.isInited()) { mcGeoActive = AnimChannel<float>(geoactive); }
189
190         checkIsAnimated();
191
192         mIsInitialized = true;
193         debMsgStd("ntlGeometryObject::initialize",DM_MSG,"GeoObj '"<<this->getName()<<"': visible="<<this->mVisible<<" gid="<<this->mGeoInitId<<" gtype="<<mGeoInitType<<","<<ginitStr<<
194                         " gvel="<<mInitialVelocity<<" gisect="<<mGeoInitIntersect, 10); // debug
195 }
196
197 /*! notify object that dump is in progress (e.g. for particles) */
198 // default action - do nothing...
199 void ntlGeometryObject::notifyOfDump(int dumtp, int frameNr,char *frameNrStr,string outfilename, double simtime) {
200   bool debugOut=false;
201         if(debugOut) debMsgStd("ntlGeometryObject::notifyOfDump",DM_MSG," dt:"<<dumtp<<" obj:"<<this->getName()<<" frame:"<<frameNrStr<<","<<frameNr<<",t"<<simtime<<" to "<<outfilename, 10); // DEBUG
202 }
203
204 /*****************************************************************************/
205 /* Search the material for this object from the material list */
206 /*****************************************************************************/
207 void ntlGeometryObject::searchMaterial(vector<ntlMaterial *> *mat)
208 {
209         /* search the list... */
210         int i=0;
211         for (vector<ntlMaterial*>::iterator iter = mat->begin();
212          iter != mat->end(); iter++) {
213                 if( mMaterialName == (*iter)->getName() ) {
214                         //warnMsg("ntlGeometryObject::searchMaterial","for obj '"<<getName()<<"' found - '"<<(*iter)->getName()<<"' "<<i); // DEBUG
215                         mpMaterial = (*iter);
216                         return;
217                 }
218                 i++;
219         }
220         errFatal("ntlGeometryObject::searchMaterial","Unknown material '"<<mMaterialName<<"' ! ", SIMWORLD_INITERROR);
221         mpMaterial = new ntlMaterial();
222         return;
223 }
224
225 /******************************************************************************
226  * static add triangle function
227  *****************************************************************************/
228 void ntlGeometryObject::sceneAddTriangle(
229                 ntlVec3Gfx  p1,ntlVec3Gfx  p2,ntlVec3Gfx  p3,
230                 ntlVec3Gfx pn1,ntlVec3Gfx pn2,ntlVec3Gfx pn3,
231                 ntlVec3Gfx trin, bool smooth,
232                 vector<ntlTriangle> *triangles,
233                 vector<ntlVec3Gfx>  *vertices,
234                 vector<ntlVec3Gfx>  *normals) {
235         ntlTriangle tri;
236         int tempVert;
237   
238         if(normals->size() != vertices->size()) {
239                 errFatal("ntlGeometryObject::sceneAddTriangle","For '"<<this->mName<<"': Vertices and normals sizes to not match!!!",SIMWORLD_GENERICERROR);
240
241         } else {
242                 
243                 vertices->push_back( p1 ); 
244                 normals->push_back( pn1 ); 
245                 tempVert = normals->size()-1;
246                 tri.getPoints()[0] = tempVert;
247                 
248                 vertices->push_back( p2 ); 
249                 normals->push_back( pn2 ); 
250                 tempVert = normals->size()-1;
251                 tri.getPoints()[1] = tempVert;
252                 
253                 vertices->push_back( p3 ); 
254                 normals->push_back( pn3 ); 
255                 tempVert = normals->size()-1;
256                 tri.getPoints()[2] = tempVert;
257                 
258                 
259                 /* init flags from ntl_ray.h */
260                 int flag = 0; 
261                 if(getVisible()){ flag |= TRI_GEOMETRY; }
262                 if(getCastShadows() ) { 
263                         flag |= TRI_CASTSHADOWS; } 
264                 
265                 /* init geo init id */
266                 int geoiId = getGeoInitId(); 
267                 //if((geoiId > 0) && (mVolumeInit&VOLUMEINIT_VOLUME) && (!mIsAnimated)) { 
268                 if((geoiId > 0) && (mVolumeInit&VOLUMEINIT_VOLUME)) { 
269                         flag |= (1<< (geoiId+4)); 
270                         flag |= mGeoInitType; 
271                 } 
272                 /*errMsg("ntlScene::addTriangle","DEBUG flag="<<convertFlags2String(flag) ); */ 
273                 tri.setFlags( flag );
274                 
275                 /* triangle normal missing */
276                 tri.setNormal( trin );
277                 tri.setSmoothNormals( smooth );
278                 tri.setObjectId( this->mObjectId );
279                 triangles->push_back( tri ); 
280         } /* normals check*/ 
281 }
282
283 void ntlGeometryObject::sceneAddTriangleNoVert(int *trips,
284                 ntlVec3Gfx trin, bool smooth,
285                 vector<ntlTriangle> *triangles) {
286         ntlTriangle tri;
287                 
288         tri.getPoints()[0] = trips[0];
289         tri.getPoints()[1] = trips[1];
290         tri.getPoints()[2] = trips[2];
291
292         // same as normal sceneAddTriangle
293
294         /* init flags from ntl_ray.h */
295         int flag = 0; 
296         if(getVisible()){ flag |= TRI_GEOMETRY; }
297         if(getCastShadows() ) { 
298                 flag |= TRI_CASTSHADOWS; } 
299
300         /* init geo init id */
301         int geoiId = getGeoInitId(); 
302         if((geoiId > 0) && (mVolumeInit&VOLUMEINIT_VOLUME)) { 
303                 flag |= (1<< (geoiId+4)); 
304                 flag |= mGeoInitType; 
305         } 
306         /*errMsg("ntlScene::addTriangle","DEBUG flag="<<convertFlags2String(flag) ); */ 
307         tri.setFlags( flag );
308
309         /* triangle normal missing */
310         tri.setNormal( trin );
311         tri.setSmoothNormals( smooth );
312         tri.setObjectId( this->mObjectId );
313         triangles->push_back( tri ); 
314 }
315
316
317 /******************************************************************************/
318 /* Init channels from float arrays (for elbeem API) */
319 /******************************************************************************/
320
321 #define ADD_CHANNEL_VEC(dst,nvals,val) \
322                 vals.clear(); time.clear(); elbeemSimplifyChannelVec3(val,&nvals); \
323                 for(int i=0; i<(nvals); i++) { \
324                         vals.push_back(ntlVec3Gfx((val)[i*4+0], (val)[i*4+1],(val)[i*4+2] )); \
325                         time.push_back( (val)[i*4+3] ); \
326                 } \
327                 (dst) = AnimChannel< ntlVec3Gfx >(vals,time); 
328
329 #define ADD_CHANNEL_FLOAT(dst,nvals,val) \
330                 valsfloat.clear(); time.clear(); elbeemSimplifyChannelFloat(val,&nvals); \
331                 for(int i=0; i<(nvals); i++) { \
332                         valsfloat.push_back( (val)[i*2+0] ); \
333                         time.push_back( (val)[i*2+1] ); \
334                 } \
335                 (dst) = AnimChannel< float >(valsfloat,time); 
336
337 void ntlGeometryObject::initChannels(
338                 int nTrans, float *trans, int nRot, float *rot, int nScale, float *scale,
339                 int nAct, float *act, int nIvel, float *ivel,
340                 int nAttrFStr, float *attrFStr,
341                 int nAttrFRad, float *attrFRad,
342                 int nVelFStr, float *velFStr,
343                 int nVelFRad, float *velFRad
344                 ) {
345         const bool debugInitc=true;
346         if(debugInitc) { debMsgStd("ntlGeometryObject::initChannels",DM_MSG,"nt:"<<nTrans<<" nr:"<<nRot<<" ns:"<<nScale, 10); 
347                          debMsgStd("ntlGeometryObject::initChannels",DM_MSG,"na:"<<nAct<<" niv:"<<nIvel<<" ", 10); }
348         vector<ntlVec3Gfx> vals;
349         vector<float>  valsfloat;
350         vector<double> time;
351         if((trans)&&(nTrans>0)) {  ADD_CHANNEL_VEC(mcTrans, nTrans, trans); }
352         if((rot)&&(nRot>0)) {      ADD_CHANNEL_VEC(mcRot, nRot, rot); }
353         if((scale)&&(nScale>0)) {  ADD_CHANNEL_VEC(mcScale, nScale, scale); }
354         if((act)&&(nAct>0)) {      ADD_CHANNEL_FLOAT(mcGeoActive, nAct, act); }
355         if((ivel)&&(nIvel>0)) {    ADD_CHANNEL_VEC(mcInitialVelocity, nIvel, ivel); }
356         
357         /* fluid control channels */
358         if((attrFStr)&&(nAttrFStr>0)) { ADD_CHANNEL_FLOAT(mcAttrFStr, nAttrFStr, attrFStr); }
359         if((attrFRad)&&(nAttrFRad>0)) { ADD_CHANNEL_FLOAT(mcAttrFRad, nAttrFRad, attrFRad); }
360         if((velFStr)&&(nVelFStr>0)) {   ADD_CHANNEL_FLOAT(mcVelFStr, nAct, velFStr); }
361         if((velFRad)&&(nVelFRad>0)) {   ADD_CHANNEL_FLOAT(mcVelFRad, nVelFRad, velFRad); }
362
363         checkIsAnimated();
364         
365         if(debugInitc) { 
366                 debMsgStd("ntlGeometryObject::initChannels",DM_MSG,getName()<<
367                                 " nt:"<<mcTrans.accessValues().size()<<" nr:"<<mcRot.accessValues().size()<<
368                                 " ns:"<<mcScale.accessValues().size()<<" isAnim:"<<mIsAnimated, 10); }
369
370         if(debugInitc) {
371                 std::ostringstream ostr;
372                 ostr << "trans: ";
373                 for(size_t i=0; i<mcTrans.accessValues().size(); i++) {
374                         ostr<<" "<<mcTrans.accessValues()[i]<<"@"<<mcTrans.accessTimes()[i]<<" ";
375                 } ostr<<";   ";
376                 ostr<<"rot: ";
377                 for(size_t i=0; i<mcRot.accessValues().size(); i++) {
378                         ostr<<" "<<mcRot.accessValues()[i]<<"@"<<mcRot.accessTimes()[i]<<" ";
379                 } ostr<<";   ";
380                 ostr<<"scale: ";
381                 for(size_t i=0; i<mcScale.accessValues().size(); i++) {
382                         ostr<<" "<<mcScale.accessValues()[i]<<"@"<<mcScale.accessTimes()[i]<<" ";
383                 } ostr<<";   ";
384                 ostr<<"act: ";
385                 for(size_t i=0; i<mcGeoActive.accessValues().size(); i++) {
386                         ostr<<" "<<mcGeoActive.accessValues()[i]<<"@"<<mcGeoActive.accessTimes()[i]<<" ";
387                 } ostr<<";   ";
388                 ostr<<"ivel: ";
389                 for(size_t i=0; i<mcInitialVelocity.accessValues().size(); i++) {
390                         ostr<<" "<<mcInitialVelocity.accessValues()[i]<<"@"<<mcInitialVelocity.accessTimes()[i]<<" ";
391                 } ostr<<";   ";
392                 debMsgStd("ntlGeometryObject::initChannels",DM_MSG,"Inited "<<ostr.str(),10);
393         }
394 }
395 #undef ADD_CHANNEL
396
397
398 /*****************************************************************************/
399 /* apply object translation at time t*/
400 /*****************************************************************************/
401 void ntlGeometryObject::applyTransformation(double t, vector<ntlVec3Gfx> *verts, vector<ntlVec3Gfx> *norms, int vstart, int vend, int forceTrafo) {
402         if(    (mcTrans.accessValues().size()>1)  // VALIDATE
403             || (mcRot.accessValues().size()>1) 
404             || (mcScale.accessValues().size()>1) 
405             || (forceTrafo)
406             || (!mHaveCachedMov)
407                 ) {
408                 // transformation is animated, continue
409                 ntlVec3Gfx pos = getTranslation(t); 
410                 ntlVec3Gfx scale = mcScale.get(t);
411                 ntlVec3Gfx rot = mcRot.get(t);
412                 ntlMat4Gfx rotMat;
413                 rotMat.initRotationXYZ(rot[0],rot[1],rot[2]);
414                 pos += mInitialPos;
415                 errMsg("ntlGeometryObject::applyTransformation","obj="<<getName()<<" t"<<pos<<" r"<<rot<<" s"<<scale);
416                 for(int i=vstart; i<vend; i++) {
417                         (*verts)[i] *= scale;
418                         (*verts)[i] = rotMat * (*verts)[i];
419                         (*verts)[i] += pos;
420                         //if(i<10) errMsg("ntlGeometryObject::applyTransformation"," v"<<i<<"/"<<vend<<"="<<(*verts)[i]);
421                 }
422                 if(norms) {
423                         for(int i=vstart; i<vend; i++) {
424                                 (*norms)[i] = rotMat * (*norms)[i];
425                         }
426                 }
427         } else {
428                 // not animated, cached points were already returned
429                 errMsg ("ntlGeometryObject::applyTransformation","Object "<<getName()<<" used cached points ");
430         }
431 }
432
433 /*! init triangle divisions */
434 void ntlGeometryObject::calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri) {
435         mTriangleDivs1.resize( tris.size() );
436         mTriangleDivs2.resize( tris.size() );
437
438         //fsTri *= 2.; // DEBUG! , wrong init!
439
440         for(size_t i=0; i<tris.size(); i++) {
441                 const ntlVec3Gfx p0 = verts[ tris[i].getPoints()[0] ];
442                 const ntlVec3Gfx p1 = verts[ tris[i].getPoints()[1] ];
443                 const ntlVec3Gfx p2 = verts[ tris[i].getPoints()[2] ];
444                 const ntlVec3Gfx side1 = p1 - p0;
445                 const ntlVec3Gfx side2 = p2 - p0;
446                 int divs1=0, divs2=0;
447                 if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
448                 if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
449
450                 mTriangleDivs1[i] = divs1;
451                 mTriangleDivs2[i] = divs2;
452         }
453 }
454
455 /*! Prepare points for moving objects */
456 void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
457         if((mMovPntsInited==featureSize)&&(!getMeshAnimated())) return;
458         const bool debugMoinit=false;
459
460         vector<ntlTriangle> triangles; 
461         vector<ntlVec3Gfx> vertices; 
462         vector<ntlVec3Gfx> vnormals; 
463         int objectId = 1;
464         this->getTriangles(time, &triangles,&vertices,&vnormals,objectId);
465         
466         mMovPoints.clear();
467         mMovNormals.clear();
468         if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" has v:"<<vertices.size()<<" t:"<<triangles.size() );
469         // no points?
470         if(vertices.size()<1) {
471                 mMaxMovPnt=-1;
472                 return; 
473         }
474         ntlVec3f maxscale = channelFindMaxVf(mcScale);
475         float maxpart = ABS(maxscale[0]);
476         if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
477         if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
478         float scaleFac = 1.0/(maxpart);
479         // TODO - better reinit from time to time?
480         const gfxReal fsTri = featureSize*0.5 *scaleFac;
481         if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","maxscale:"<<maxpart<<" featureSize:"<<featureSize<<" fsTri:"<<fsTri );
482
483         if(mTriangleDivs1.size()!=triangles.size()) {
484                 calcTriangleDivs(vertices,triangles,fsTri);
485         }
486
487         // debug: count points to init
488         /*if(debugMoinit) {
489                 errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" estimating...");
490                 int countp=vertices.size()*2;
491                 for(size_t i=0; i<triangles.size(); i++) {
492                         ntlVec3Gfx p0 = vertices[ triangles[i].getPoints()[0] ];
493                         ntlVec3Gfx side1 = vertices[ triangles[i].getPoints()[1] ] - p0;
494                         ntlVec3Gfx side2 = vertices[ triangles[i].getPoints()[2] ] - p0;
495                         int divs1=0, divs2=0;
496                         if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
497                         if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
498                         errMsg("ntlGeometryObject::initMovingPoints","tri:"<<i<<" p:"<<p0<<" s1:"<<side1<<" s2:"<<side2<<" -> "<<divs1<<","<<divs2 );
499                         if(divs1+divs2 > 0) {
500                                 for(int u=0; u<=divs1; u++) {
501                                         for(int v=0; v<=divs2; v++) {
502                                                 const gfxReal uf = (gfxReal)(u+TRI_UVOFFSET) / (gfxReal)(divs1+0.0);
503                                                 const gfxReal vf = (gfxReal)(v+TRI_UVOFFSET) / (gfxReal)(divs2+0.0);
504                                                 if(uf+vf>1.0) continue;
505                                                 countp+=2;
506                                         }
507                                 }
508                         }
509                 }
510                 errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" requires:"<<countp*2);
511         } // */
512
513         bool discardInflowBack = false;
514         if( (mGeoInitType==FGI_MBNDINFLOW) && (mcInitialVelocity.accessValues().size()<1) ) discardInflowBack = true;
515         discardInflowBack = false; // DEBUG disable for now
516
517
518         // init std points
519         for(size_t i=0; i<vertices.size(); i++) {
520                 ntlVec3Gfx p = vertices[ i ];
521                 ntlVec3Gfx n = vnormals[ i ];
522                 // discard inflow backsides
523                 //if( (mGeoInitType==FGI_MBNDINFLOW) && (!mIsAnimated)) {
524                 if(discardInflowBack) { //if( (mGeoInitType==FGI_MBNDINFLOW) && (!mIsAnimated)) {
525                         if(dot(mInitialVelocity,n)<0.0) continue;
526                 }
527                 mMovPoints.push_back(p);
528                 mMovNormals.push_back(n);
529                 if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","std"<<i<<" p"<<p<<" n"<<n<<" ");
530         }
531         // init points & refine...
532         for(size_t i=0; i<triangles.size(); i++) {
533                 int *trips = triangles[i].getPoints();
534                 const ntlVec3Gfx p0 = vertices[ trips[0] ];
535                 const ntlVec3Gfx side1 = vertices[ trips[1] ] - p0;
536                 const ntlVec3Gfx side2 = vertices[ trips[2] ] - p0;
537                 int divs1=mTriangleDivs1[i], divs2=mTriangleDivs2[i];
538                 
539                 const ntlVec3Gfx trinormOrg = getNormalized(cross(side1,side2));
540                 const ntlVec3Gfx trinorm = trinormOrg*0.25*featureSize;
541                 if(discardInflowBack) { 
542                         if(dot(mInitialVelocity,trinorm)<0.0) continue;
543                 }
544                 if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","Tri1 "<<vertices[trips[0]]<<","<<vertices[trips[1]]<<","<<vertices[trips[2]]<<" "<<divs1<<","<<divs2 );
545                 if(divs1+divs2 > 0) {
546                         for(int u=0; u<=divs1; u++) {
547                                 for(int v=0; v<=divs2; v++) {
548                                         const gfxReal uf = (gfxReal)(u+TRI_UVOFFSET) / (gfxReal)(divs1+0.0);
549                                         const gfxReal vf = (gfxReal)(v+TRI_UVOFFSET) / (gfxReal)(divs2+0.0);
550                                         if(uf+vf>1.0) continue;
551                                         ntlVec3Gfx p = 
552                                                 vertices[ trips[0] ] * (1.0-uf-vf)+
553                                                 vertices[ trips[1] ] * uf +
554                                                 vertices[ trips[2] ] * vf;
555                                         //ntlVec3Gfx n = vnormals[ 
556                                                 //trips[0] ] * (1.0-uf-vf)+
557                                                 //vnormals[ trips[1] ]*uf +
558                                                 //vnormals[ trips[2] ]*vf;
559                                         //normalize(n);
560                                         // discard inflow backsides
561
562                                         mMovPoints.push_back(p + trinorm); // NEW!?
563                                         mMovPoints.push_back(p - trinorm);
564                                         mMovNormals.push_back(trinormOrg);
565                                         mMovNormals.push_back(trinormOrg); 
566                                         //errMsg("TRINORM","p"<<p<<" n"<<n<<" trin"<<trinorm);
567                                 }
568                         }
569                 }
570         }
571
572         // find max point
573         mMaxMovPnt = 0;
574         gfxReal dist = normNoSqrt(mMovPoints[0]);
575         for(size_t i=0; i<mMovPoints.size(); i++) {
576                 if(normNoSqrt(mMovPoints[i])>dist) {
577                         mMaxMovPnt = i;
578                         dist = normNoSqrt(mMovPoints[0]);
579                 }
580         }
581
582         if(    (this->getMeshAnimated())
583       || (mcTrans.accessValues().size()>1)  // VALIDATE
584             || (mcRot.accessValues().size()>1) 
585             || (mcScale.accessValues().size()>1) 
586                 ) {
587                 // also do trafo...
588         } else {
589                 mCachedMovPoints = mMovPoints;
590                 mCachedMovNormals = mMovNormals;
591                 //applyTransformation(time, &mCachedMovPoints, &mCachedMovNormals, 0, mCachedMovPoints.size(), true);
592                 applyTransformation(time, &mCachedMovPoints, NULL, 0, mCachedMovPoints.size(), true);
593                 mHaveCachedMov = true;
594                 debMsgStd("ntlGeometryObject::initMovingPoints",DM_MSG,"Object "<<getName()<<" cached points ", 7);
595         }
596
597         mMovPntsInited = featureSize;
598         debMsgStd("ntlGeometryObject::initMovingPoints",DM_MSG,"Object "<<getName()<<" inited v:"<<vertices.size()<<"->"<<mMovPoints.size() , 5);
599 }
600 /*! Prepare points for animated objects,
601  * init both sets, never used cached points 
602  * discardInflowBack ignore */
603 void ntlGeometryObject::initMovingPointsAnim(
604                  double srctime, vector<ntlVec3Gfx> &srcmovPoints,
605                  double dsttime, vector<ntlVec3Gfx> &dstmovPoints,
606                  vector<ntlVec3Gfx> *dstmovNormals,
607                  gfxReal featureSize,
608                  ntlVec3Gfx geostart, ntlVec3Gfx geoend
609                  ) {
610         const bool debugMoinit=false;
611
612         vector<ntlTriangle> srctriangles; 
613         vector<ntlVec3Gfx> srcvertices; 
614         vector<ntlVec3Gfx> unused_normals; 
615         vector<ntlTriangle> dsttriangles; 
616         vector<ntlVec3Gfx> dstvertices; 
617         vector<ntlVec3Gfx> dstnormals; 
618         int objectId = 1;
619         // TODO optimize? , get rid of normals?
620         unused_normals.clear();
621         this->getTriangles(srctime, &srctriangles,&srcvertices,&unused_normals,objectId);
622         unused_normals.clear();
623         this->getTriangles(dsttime, &dsttriangles,&dstvertices,&dstnormals,objectId);
624         
625         srcmovPoints.clear();
626         dstmovPoints.clear();
627         if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has srcv:"<<srcvertices.size()<<" srct:"<<srctriangles.size() );
628         if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has dstv:"<<dstvertices.size()<<" dstt:"<<dsttriangles.size() );
629         // no points?
630         if(srcvertices.size()<1) {
631                 mMaxMovPnt=-1;
632                 return; 
633         }
634         if((srctriangles.size() != dsttriangles.size()) ||
635            (srcvertices.size() != dstvertices.size()) ) {
636                 errMsg("ntlGeometryObject::initMovingPointsAnim","Invalid triangle numbers! Aborting...");
637                 return;
638         }
639         ntlVec3f maxscale = channelFindMaxVf(mcScale);
640         float maxpart = ABS(maxscale[0]);
641         if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
642         if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
643         float scaleFac = 1.0/(maxpart);
644         // TODO - better reinit from time to time?
645         const gfxReal fsTri = featureSize*0.5 *scaleFac;
646         if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","maxscale:"<<maxpart<<" featureSize:"<<featureSize<<" fsTri:"<<fsTri );
647
648         if(mTriangleDivs1.size()!=srctriangles.size()) {
649                 calcTriangleDivs(srcvertices,srctriangles,fsTri);
650         }
651
652
653         // init std points
654         for(size_t i=0; i<srcvertices.size(); i++) {
655                 srcmovPoints.push_back(srcvertices[i]);
656                 //srcmovNormals.push_back(srcnormals[i]);
657         }
658         for(size_t i=0; i<dstvertices.size(); i++) {
659                 dstmovPoints.push_back(dstvertices[i]);
660                 if(dstmovNormals) (*dstmovNormals).push_back(dstnormals[i]);
661         }
662         if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<" dst:"<<dstmovPoints.size()<<" " );
663         // init points & refine...
664         for(size_t i=0; i<srctriangles.size(); i++) {
665                 const int divs1=mTriangleDivs1[i];
666                 const int       divs2=mTriangleDivs2[i];
667                 if(divs1+divs2 > 0) {
668                         int *srctrips = srctriangles[i].getPoints();
669                         int *dsttrips = dsttriangles[i].getPoints();
670                         const ntlVec3Gfx srcp0 =    srcvertices[ srctrips[0] ];
671                         const ntlVec3Gfx srcside1 = srcvertices[ srctrips[1] ] - srcp0;
672                         const ntlVec3Gfx srcside2 = srcvertices[ srctrips[2] ] - srcp0;
673                         const ntlVec3Gfx dstp0 =    dstvertices[ dsttrips[0] ];
674                         const ntlVec3Gfx dstside1 = dstvertices[ dsttrips[1] ] - dstp0;
675                         const ntlVec3Gfx dstside2 = dstvertices[ dsttrips[2] ] - dstp0;
676                         const ntlVec3Gfx src_trinorm    = getNormalized(cross(srcside1,srcside2))*0.25*featureSize;
677                         const ntlVec3Gfx dst_trinormOrg = getNormalized(cross(dstside1,dstside2));
678                         const ntlVec3Gfx dst_trinorm    = dst_trinormOrg                         *0.25*featureSize;
679                         //errMsg("ntlGeometryObject::initMovingPointsAnim","Tri1 "<<srcvertices[srctrips[0]]<<","<<srcvertices[srctrips[1]]<<","<<srcvertices[srctrips[2]]<<" "<<divs1<<","<<divs2 );
680                         for(int u=0; u<=divs1; u++) {
681                                 for(int v=0; v<=divs2; v++) {
682                                         const gfxReal uf = (gfxReal)(u+TRI_UVOFFSET) / (gfxReal)(divs1+0.0);
683                                         const gfxReal vf = (gfxReal)(v+TRI_UVOFFSET) / (gfxReal)(divs2+0.0);
684                                         if(uf+vf>1.0) continue;
685                                         ntlVec3Gfx srcp = 
686                                                 srcvertices[ srctrips[0] ] * (1.0-uf-vf)+
687                                                 srcvertices[ srctrips[1] ] * uf +
688                                                 srcvertices[ srctrips[2] ] * vf;
689                                         ntlVec3Gfx dstp = 
690                                                 dstvertices[ dsttrips[0] ] * (1.0-uf-vf)+
691                                                 dstvertices[ dsttrips[1] ] * uf +
692                                                 dstvertices[ dsttrips[2] ] * vf;
693
694                                         // cutoffDomain
695                                         if((srcp[0]<geostart[0]) && (dstp[0]<geostart[0])) continue;
696                                         if((srcp[1]<geostart[1]) && (dstp[1]<geostart[1])) continue;
697                                         if((srcp[2]<geostart[2]) && (dstp[2]<geostart[2])) continue;
698                                         if((srcp[0]>geoend[0]  ) && (dstp[0]>geoend[0]  )) continue;
699                                         if((srcp[1]>geoend[1]  ) && (dstp[1]>geoend[1]  )) continue;
700                                         if((srcp[2]>geoend[2]  ) && (dstp[2]>geoend[2]  )) continue;
701                                         
702                                         srcmovPoints.push_back(srcp+src_trinorm); // SURFENHTEST
703                                         srcmovPoints.push_back(srcp-src_trinorm);
704
705                                         dstmovPoints.push_back(dstp+dst_trinorm); // SURFENHTEST
706                                         dstmovPoints.push_back(dstp-dst_trinorm);
707                                         if(dstmovNormals) {
708                                                 (*dstmovNormals).push_back(dst_trinormOrg);
709                                                 (*dstmovNormals).push_back(dst_trinormOrg); }
710                                 }
711                         }
712                 }
713         }
714
715         // find max point not necessary
716         debMsgStd("ntlGeometryObject::initMovingPointsAnim",DM_MSG,"Object "<<getName()<<" inited v:"<<srcvertices.size()<<"->"<<srcmovPoints.size()<<","<<dstmovPoints.size() , 5);
717 }
718
719 /*! Prepare points for moving objects */
720 void ntlGeometryObject::getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3Gfx> *norms) {
721         if(mHaveCachedMov) {
722                 ret = mCachedMovPoints;
723                 if(norms) { 
724                         *norms = mCachedMovNormals; 
725                         //errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
726                 }
727                 //errMsg ("ntlGeometryObject::getMovingPoints","Object "<<getName()<<" used cached points "); // DEBUG
728                 return;
729         }
730
731         ret = mMovPoints;
732         if(norms) { 
733                 //errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
734                 *norms = mMovNormals; 
735         }
736 }
737
738
739 /*! Calculate max. velocity on object from t1 to t2 */
740 ntlVec3Gfx ntlGeometryObject::calculateMaxVel(double t1, double t2) {
741         ntlVec3Gfx vel(0.);
742         if(mMaxMovPnt<0) return vel;
743                 
744         vector<ntlVec3Gfx> verts1,verts2;
745         verts1.push_back(mMovPoints[mMaxMovPnt]);
746         verts2 = verts1;
747         applyTransformation(t1,&verts1,NULL, 0,verts1.size(), true);
748         applyTransformation(t2,&verts2,NULL, 0,verts2.size(), true);
749
750         vel = (verts2[0]-verts1[0]); // /(t2-t1);
751         //errMsg("ntlGeometryObject::calculateMaxVel","t1="<<t1<<" t2="<<t2<<" p1="<<verts1[0]<<" p2="<<verts2[0]<<" v="<<vel);
752         return vel;
753 }
754
755 /*! get translation at time t*/
756 ntlVec3Gfx ntlGeometryObject::getTranslation(double t) {
757         ntlVec3Gfx pos = mcTrans.get(t);
758   // DEBUG CP_FORCECIRCLEINIT 1
759         /*if( 
760                         (mName.compare(string("0__ts1"))==0) ||
761                         (mName.compare(string("1__ts1"))==0) ||
762                         (mName.compare(string("2__ts1"))==0) ||
763                         (mName.compare(string("3__ts1"))==0) ||
764                         (mName.compare(string("4__ts1"))==0) ||
765                         (mName.compare(string("5__ts1"))==0) ||
766                         (mName.compare(string("6__ts1"))==0) ||
767                         (mName.compare(string("7__ts1"))==0) ||
768                         (mName.compare(string("8__ts1"))==0) ||
769                         (mName.compare(string("9__ts1"))==0) 
770                         ) { int j=mName[0]-'0';
771                         ntlVec3Gfx ppos(0.); { // DEBUG
772                         const float tscale=10.;
773                         const float tprevo = 0.33;
774                         const ntlVec3Gfx   toff(50,50,0);
775                         const ntlVec3Gfx oscale(30,30,0);
776                         ppos[0] =  cos(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[0] + toff[0];
777                         ppos[1] = -sin(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[1] + toff[1];
778                         ppos[2] =                               toff[2]; } // DEBUG
779                         pos = ppos;
780                         pos[2] = 0.15;
781         }
782   // DEBUG CP_FORCECIRCLEINIT 1 */
783         return pos;
784 }
785 /*! get active flag time t*/
786 float ntlGeometryObject::getGeoActive(double t) {
787         float act = (mcGeoActive.get(t) >= 1.) ? 1.0 : 0.0; 
788         return act;
789 }
790
791 void ntlGeometryObject::setInitialVelocity(ntlVec3Gfx set) { 
792         mInitialVelocity=set; 
793         mcInitialVelocity = AnimChannel<ntlVec3Gfx>(set); 
794 }
795 ntlVec3Gfx ntlGeometryObject::getInitialVelocity(double t) { 
796         ntlVec3Gfx v =  mcInitialVelocity.get(t); //return mInitialVelocity; 
797         if(!mLocalCoordInivel) return v;
798
799         ntlVec3Gfx rot = mcRot.get(t);
800         ntlMat4Gfx rotMat;
801         rotMat.initRotationXYZ(rot[0],rot[1],rot[2]);
802         v = rotMat * v;
803         return v;
804 }
805         
806