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