This version now includes the fluid control sources, however the Blender
[blender.git] / intern / elbeem / intern / mvmcoords.cpp
1 /******************************************************************************
2  *
3 // El'Beem - the visual lattice boltzmann freesurface simulator
4 // All code distributed as part of El'Beem is covered by the version 2 of the 
5 // GNU General Public License. See the file COPYING for details.  
6 //
7 // Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
8 //
9  *
10  * Mean Value Mesh Coords class
11  *
12  *****************************************************************************/
13
14 #include "mvmcoords.h"
15 #include <algorithm>
16 using std::vector;
17
18 void MeanValueMeshCoords::clear() 
19 {
20         mVertices.resize(0);
21         mNumVerts = 0;
22 }
23
24 void MeanValueMeshCoords::calculateMVMCs(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle> &tris, 
25                 vector<ntlVec3Gfx> &points, gfxReal numweights)
26 {
27         clear();
28         mvmTransferPoint tds;
29         int mem = 0;
30         int i = 0;
31
32         mNumVerts = (int)reference_vertices.size();
33
34         for (vector<ntlVec3Gfx>::iterator iter = points.begin(); iter != points.end(); ++iter, ++i) {
35                 if(i%(points.size()/10)==1) debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"Computing weights, points: "<<i<<"/"<<points.size(),5 );
36                 tds.lastpos = *iter;
37                 tds.weights.resize(0);  // clear
38                 computeWeights(reference_vertices, tris, tds, numweights);
39                 mem += (int)tds.weights.size();
40                 mVertices.push_back(tds);
41         }    
42         int mbmem = mem * sizeof(mvmFloat) / (1024*1024);
43         debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"vertices:"<<mNumVerts<<" points:"<<points.size()<<" weights:"<<mem<<", wmem:"<<mbmem<<"MB ",7 );
44 }
45
46 // from: mean value coordinates for closed triangular meshes
47 // attention: fails if a point is exactly (or very close) to a vertex
48 void MeanValueMeshCoords::computeWeights(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle>& tris,
49                 mvmTransferPoint& tds, gfxReal numweights)
50 {
51         const bool mvmFullDebug=false;
52         //const ntlVec3Gfx cEPS = 1.0e-6;
53         const mvmFloat cEPS = 1.0e-14;
54
55         //mvmFloat d[3], s[3], phi[3],c[3];
56         ntlVec3d u[3],c,d,s,phi;
57         int indices[3];
58
59         for (int i = 0; i < (int)reference_vertices.size(); ++i) {
60                 tds.weights.push_back(mvmIndexWeight(i, 0.0));
61         }
62
63         // for each triangle
64         //for (vector<ntlTriangle>::iterator iter = tris.begin(); iter != tris.end();) {
65         for(int t=0; t<(int)tris.size(); t++) {
66
67                 for (int i = 0; i < 3; ++i) { //, ++iter) {
68                         indices[i] = tris[t].getPoints()[i];
69                         u[i] = vec2D(reference_vertices[ indices[i] ]-tds.lastpos);
70                         d[i] = normalize(u[i]); //.normalize();        
71                         //assert(d[i] != 0.);
72                         if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","t"<<t<<" i"<<indices[i] //<<" lp"<<tds.lastpos
73                                         <<" v"<<reference_vertices[indices[i]]<<" u"<<u[i]<<" ");
74                         // on vertex!
75                         //? if(d[i]<=0.) continue;
76                 }
77                 //for (int i = 0; i < 3; ++i) { errMsg("III"," "<<i             <<" i"<<indices[i]<<reference_vertices[ indices[i] ] ); }
78
79                 // arcsin is not needed, see paper
80                 phi[0] = 2.*asin( (mvmFloat)(0.5* norm(u[1]-u[2]) )  );
81                 phi[1] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[2]) )  );
82                 phi[2] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[1]) )  );
83                 mvmFloat h = (phi[0] + phi[1] + phi[2])*0.5;
84                 if (M_PI-h < cEPS) {
85                         if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","point on triangle");
86                         tds.weights.resize(0);
87                         tds.weights.push_back( mvmIndexWeight(indices[0], sin(phi[0])*d[1]*d[2]));
88                         tds.weights.push_back( mvmIndexWeight(indices[1], sin(phi[1])*d[0]*d[2]));
89                         tds.weights.push_back( mvmIndexWeight(indices[2], sin(phi[2])*d[1]*d[0]));
90                         break;
91                 }
92                 mvmFloat sinh = 2.*sin(h);
93                 c[0] = (sinh*sin(h-phi[0]))/(sin(phi[1])*sin(phi[2]))-1.;
94                 c[1] = (sinh*sin(h-phi[1]))/(sin(phi[0])*sin(phi[2]))-1.;    
95                 c[2] = (sinh*sin(h-phi[2]))/(sin(phi[0])*sin(phi[1]))-1.;   
96                 if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","c="<<c<<" phi="<<phi<<" d="<<d);
97                 //if (c[0] > 1. || c[0] < 0. || c[1] > 1. || c[1] < 0. || c[2] > 1. || c[2] < 0.) continue;
98
99                 s[0] = sqrtf((float)(1.-c[0]*c[0]));
100                 s[1] = sqrtf((float)(1.-c[1]*c[1]));
101                 s[2] = sqrtf((float)(1.-c[2]*c[2]));
102
103                 if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","s");
104                 if (s[0] <= cEPS || s[1] <= cEPS || s[2] <= cEPS) {
105                         //MSG("position lies outside the triangle on the same plane -> ignore it");
106                         continue;
107                 }
108                 const mvmFloat u0x = u[0][0];
109                 const mvmFloat u0y = u[0][1];
110                 const mvmFloat u0z = u[0][2];
111                 const mvmFloat u1x = u[1][0];
112                 const mvmFloat u1y = u[1][1];
113                 const mvmFloat u1z = u[1][2];
114                 const mvmFloat u2x = u[2][0];
115                 const mvmFloat u2y = u[2][1];
116                 const mvmFloat u2z = u[2][2];
117                 mvmFloat det = u0x*u1y*u2z - u0x*u1z*u2y + u0y*u1z*u2x - u0y*u1x*u2z + u0z*u1x*u2y - u0z*u1y*u2x;
118                 //assert(det != 0.);
119                 if (det < 0.) {
120                         s[0] = -s[0];
121                         s[1] = -s[1];
122                         s[2] = -s[2];
123                 }
124
125                 tds.weights[indices[0]].weight += (phi[0]-c[1]*phi[2]-c[2]*phi[1])/(d[0]*sin(phi[1])*s[2]);
126                 tds.weights[indices[1]].weight += (phi[1]-c[2]*phi[0]-c[0]*phi[2])/(d[1]*sin(phi[2])*s[0]);
127                 tds.weights[indices[2]].weight += (phi[2]-c[0]*phi[1]-c[1]*phi[0])/(d[2]*sin(phi[0])*s[1]);
128                 if(mvmFullDebug) { errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[0]<<" o"<<tds.weights[indices[0]].weight);
129                 errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[1]<<" o"<<tds.weights[indices[1]].weight);
130                 errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[2]<<" o"<<tds.weights[indices[2]].weight);
131                 errMsg("MeanValueMeshCoords::computeWeights","\n\n\n"); }
132         }
133
134         //sort weights
135         if((numweights>0.)&& (numweights<1.) ) {
136         //if( ((int)tds.weights.size() > maxNumWeights) && (maxNumWeights > 0) ) {
137           int maxNumWeights = (int)(tds.weights.size()*numweights);
138                 if(maxNumWeights<=0) maxNumWeights = 1;
139                 std::sort(tds.weights.begin(), tds.weights.end(), std::greater<mvmIndexWeight>());
140                 // only use maxNumWeights-th largest weights
141                 tds.weights.resize(maxNumWeights);
142         }
143
144         // normalize weights
145         mvmFloat totalWeight = 0.;
146         for (vector<mvmIndexWeight>::const_iterator witer = tds.weights.begin();
147                         witer != tds.weights.end(); ++witer) {
148                 totalWeight += witer->weight;
149         }
150         mvmFloat invTotalWeight;
151         if (totalWeight == 0.) {
152                 if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","totalWeight == 0");
153                 invTotalWeight = 0.0;
154         } else {
155                 invTotalWeight = 1.0/totalWeight;
156         }
157
158         for (vector<mvmIndexWeight>::iterator viter = tds.weights.begin();
159                         viter != tds.weights.end(); ++viter) {
160                 viter->weight *= invTotalWeight;  
161                 //assert(finite(viter->weight) != 0);
162                 if(!finite(viter->weight)) viter->weight=0.;
163         }
164 }
165
166 void MeanValueMeshCoords::transfer(vector<ntlVec3Gfx> &vertices, vector<ntlVec3Gfx>& displacements) 
167 {
168         displacements.resize(0);
169
170         //debMsgStd("MeanValueMeshCoords::transfer",DM_MSG,"vertices:"<<mNumVerts<<" curr_verts:"<<vertices.size()<<" ",7 );
171         if((int)vertices.size() != mNumVerts) {
172                 errMsg("MeanValueMeshCoords::transfer","Different no of verts: "<<vertices.size()<<" vs "<<mNumVerts);
173                 return;
174         }
175
176         for (vector<mvmTransferPoint>::iterator titer = mVertices.begin(); titer != mVertices.end(); ++titer) {
177                 mvmTransferPoint &tds = *titer;
178                 ntlVec3Gfx newpos(0.0);
179
180                 for (vector<mvmIndexWeight>::iterator witer = tds.weights.begin();
181                                 witer != tds.weights.end(); ++witer) {
182                         newpos += vertices[witer->index] * witer->weight;
183                         //errMsg("transfer","np"<<newpos<<" v"<<vertices[witer->index]<<" w"<< witer->weight);
184                 }
185
186                 displacements.push_back(newpos);
187                 //displacements.push_back(newpos - tds.lastpos);
188                 //tds.lastpos = newpos;
189         }
190 }
191