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