1# (C)2007-2008 by Marc Roeder, 2# distribute under the terms of the GPL version 2.0 or later 3 4 5############################################################################# 6## This file provides functions to generate fundamental domains for 7## 3 dimensional Bieberbach groups and view them as colored polytopes in 8## JavaView. 9## 10## The functions you might want to use are in 3dimBieberbachFD.gap 11## Look there first. 12############################################################################# 13 14 15############################################################################# 16## 17## A square root approximation using continued fractions. 18## Used for the Cholesky decomposition. 19 20 21 sqrt:=function(rat) 22 local x,denom, num; 23 denom:=DenominatorRat(rat); 24 num:=NumeratorRat(rat); 25 x:=Indeterminate(Rationals); 26 return ContinuedFractionApproximationOfRoot(denom*x^2-num,100); 27 end; 28 29 30############################################################################# 31## 32## A quick and dirty implementation of a Cholesky decomposition (approximation) 33## 34 35CholeskyDecomposition:=function(mat) 36 local dim, L, line,x, col; 37 38 if not mat=TransposedMat(mat) 39 then 40 Error("matrix is not symmetric"); 41 fi; 42 dim:=DimensionSquareMat(mat); 43 if not ForAll([1..dim],d->Determinant(mat{[1..d]}{[1..d]})>0) 44 then 45 Error("matrix is not positive definite"); 46 fi; 47 L:=NullMat(dim,dim); 48 L[1][1]:=sqrt(mat[1][1]); 49 for line in [1..dim] 50 do 51 for col in [1..line-1] 52 do 53 if col>1 54 then 55 L[line][col]:=(mat[line][col]-(L[line]{[1..col-1]}*L[col]{[1..col-1]}))/L[col][col]; 56 else 57 L[line][col]:=mat[line][col]/L[col][col]; 58 fi; 59 if line>1 60 then 61 L[line][line]:=sqrt(mat[line][line]-L[line]{[1..line-1]}^2); 62 fi; 63 od; 64 od; 65 return L; 66end; 67 68 69 70############################################################################# 71## convert hsv color space to rgb color space. 72## hsv is a triple [h,s,v] where 73## h\in [0..259], s,v\in[0..255]. 74## 75hsv2rgb:=function(hsv) 76 local H, S, V, Hi, f, p, q, t, rgb; 77 78 H:=hsv[1]; 79 S:=1/256*hsv[2]; 80 V:=1/256*hsv[3]; 81 Hi:=Int(H/60) mod 6; 82 f:=H/60-Hi; 83 p:=V*(1-S); 84 q:=V*(1-f*S); 85 t:=V*(1-(1-f)*S); 86 if Hi=0 87 then 88 rgb:=[V,t,p]; 89 elif Hi=1 90 then 91 rgb:= [q,V,p]; 92 elif Hi=2 93 then 94 rgb:= [p,V,t]; 95 elif Hi=3 96 then 97 rgb:= [p,q,V]; 98 elif Hi=4 99 then 100 rgb:= [t,p,V]; 101 elif Hi=5 102 then 103 rgb:= [V,p,q]; 104 else 105 Error("bad programming, Marc"); 106 fi; 107 return List(rgb*256,Int); 108end; 109 110############################################################################# 111## create a list of length n which contains strings of the form 112## "xxx yyy zzz" representing a point in rgb-colorspace. 113## 114ncolorStrings:=function(n) 115 local colorlist, Hvalues, color, rem, H, S, V; 116 117 if n<=10 118 then 119 colorlist:=List([0..n-1]*360/n,i->[Int(i),255,255]); 120 else 121 colorlist:=[]; 122 Hvalues:=Reversed([0..Int(n/3)+1]*360/(Int(n/3)+1)); 123 for color in [1..n] 124 do 125# H:=(360/(Int(n/4)+1))*Int(color/4); 126 rem:=(color-1) mod 3; 127 if rem=0 128 then 129 H:=Remove(Hvalues); 130 S:=255; 131 V:=255; 132 elif rem=1 133 then 134 S:=127; 135 V:=255; 136 else 137 S:=255; 138 V:=127; 139 fi; 140 Add(colorlist,[H,S,V]); 141 od; 142 fi; 143 Apply(colorlist,hsv2rgb); 144 Apply(colorlist,c->List(c,String)); 145 Apply(colorlist,c->JoinStringsWithSeparator(c," ")); 146 return colorlist; 147end; 148 149 150########################## 151## the same, returning a slightly paler set of colors. 152## 153ncolorStringsPale:=function(n) 154 local colorlist, Hvalues, color, rem, H, S, V; 155 156 if n<=10 157 then 158 colorlist:=List([0..n-1]*360/n,i->[Int(i),180,200]); 159 else 160 colorlist:=[]; 161 Hvalues:=Reversed([0..Int(n/3)+1]*360/(Int(n/3)+1)); 162 for color in [1..n] 163 do 164 # H:=(360/(Int(n/4)+1))*Int(color/4); 165 rem:=(color-1) mod 3; 166 if rem=0 167 then 168 H:=Remove(Hvalues); 169 S:=127; 170 V:=255; 171 elif rem=1 172 then 173 S:=64; 174 V:=255; 175 else 176 S:=127; 177 V:=127; 178 fi; 179 Add(colorlist,[H,S,V]); 180 od; 181 fi; 182 Apply(colorlist,hsv2rgb); 183 Apply(colorlist,c->List(c,String)); 184 Apply(colorlist,c->JoinStringsWithSeparator(c," ")); 185 return colorlist; 186end; 187 188 189 190############################################################################# 191## Convert a rational number into a "floating point" string. 192## <precision> gives the number of digits after the first non-zero 193## digit after the point. 194## Example with precision=3: 1/7 is converted to "0.142" 195## 1/700 is converted to "0.00142" 196## 197 198fraction2floatString:=function(q,precision) 199 local sign, signString, magnitude, beforepoint, qright, 200 prezeros, afterpoint, returnstring; 201 202 if q=0 203 then 204 return "0"; 205 fi; 206 sign:=SignRat(q); 207 if sign=-1 208 then 209 signString:="-"; 210 else 211 signString:=""; 212 fi; 213 214 if AbsoluteValue(q)>0 215 then 216 beforepoint:=String(sign*Int(q)); 217 qright:=sign*(q-Int(q)); 218 else 219 qright:=sign*q; 220 beforepoint:="0"; 221 fi; 222 223 if qright<>0 224 then 225 magnitude:=LogInt(NumeratorRat(qright),10)-LogInt(DenominatorRat(qright),10); 226 if AbsoluteValue(qright)*10^(-magnitude)<1 227 then 228 magnitude:=magnitude-1; 229 fi; 230 prezeros:=Concatenation(List([1..-magnitude-1],i->"0")); 231 afterpoint:=String(Int(10^(precision-magnitude)*qright)); 232 returnstring:=Concatenation([signString,beforepoint,".",prezeros,afterpoint]); 233 else 234 returnstring:=Concatenation(signString,beforepoint,".0"); 235 fi; 236 while returnstring[Size(returnstring)]='0' 237 do 238 Unbind(returnstring[Size(returnstring)]); 239 od; 240 if returnstring[Size(returnstring)]='.' 241 then 242 Unbind(returnstring[Size(returnstring)]); 243 fi; 244 return returnstring; 245end; 246 247 248############################################################################# 249## This converts a number to a string and fills it up with leading zeros 250## if necessary. 251## Example with digits=3: 5->"005" 252## 1234->"1234" 253## 254numberWithLeadingZeros:=function(n,digits) 255 local string; 256 257 string:=String(n); 258 while Size(string)<digits 259 do 260 string:=Concatenation(["0",string]); 261 od; 262 return string; 263end; 264 265 266 267################################### 268 269vertexOrbitDecomposition:=function(vertexlist,group) 270 local vertices, vertexorbits, vertex, orbit; 271 272 vertices:=Set(vertexlist); 273 vertexorbits:=[]; 274 while vertices<>[] 275 do 276 vertex:=vertices[1]; 277 orbit:=Concatenation( 278 OrbitPartInVertexSetsStandardSpaceGroup(group,[vertex], 279 vertices)); 280 Add(vertexorbits,orbit); 281 SubtractSet(vertices,orbit); 282 od; 283 Apply(vertexorbits,i->List(i,j->Position(vertexlist,j))); 284 285 return vertexorbits; 286end; 287 288################################### 289 290edgeOrbitDecomposition:=function(edges,vertexlist,group) 291 local edgesasvectors, edgeorbits, edge, orbit; 292 293 edgesasvectors:=Set(edges,i->Set(i,j->vertexlist[j])); 294 edgeorbits:=[]; 295 while edgesasvectors<>[] 296 do 297 edge:=edgesasvectors[1]; 298 orbit:=Set(OrbitPartInVertexSetsStandardSpaceGroup(group,edge, 299 Set(vertexlist))); 300 SubtractSet(edgesasvectors,orbit); 301 Add(edgeorbits,orbit); 302 od; 303 edgeorbits:=Set(edgeorbits); 304 Apply(edgeorbits,i->Set(i,j->Set(j,k->Position(vertexlist,k)))); 305 return edgeorbits; 306end; 307 308 309############################################################################# 310## Take a string and put <tag> and </tag> arund it 311## 312enclosedByTag:=function(string,tag) 313 return Concatenation(["<",tag,">",string,"</",tag,">\n"]); 314end; 315 316enclosedByTagNN:=function(string,tag) 317 return Concatenation(["<",tag,">",string,"</",tag,">"]); 318end; 319 320 321enclosedByTagWithOptions:=function(string,tag,options) 322 return Concatenation(["<",tag," ",options,">\n",string,"</",tag,">\n"]); 323end; 324 325#################### 326 327javaviewWrappedDatastring:=function(title,abstract,detail,datastring,wrapper) 328 local outstring, stream, line; 329 330 outstring:=[]; 331 stream:=InputTextFile(wrapper); 332 while not IsEndOfStream(stream) 333 do 334 line:=ReadLine(stream); 335 if line="##Insert Title##\n" 336 then 337 line:=enclosedByTag(title,"title"); 338 Append(line,"\n"); 339 elif line ="##Insert Abstract##\n" 340 then 341 line:=enclosedByTag(abstract,"abstract"); 342 Append(line,"\n"); 343 344 elif line ="##Insert Detail##\n" 345 then 346 line:=enclosedByTag(detail,"detail"); 347 Append(line,"\n"); 348 elif line="##Insert Data String##\n" 349 then 350 line:=enclosedByTag(datastring,"geometries"); 351 Append(line,"\n"); 352 else 353 fi; 354 Append(outstring,line); 355 od; 356 CloseStream(stream); 357 return outstring; 358end; 359 360#################### 361 362 isposdef:=function(mat) 363 return ForAll([1..Size(mat)],d->Determinant(mat{[1..d]}{[1..d]})>0); 364 end; 365 366 367############################################################################# 368## This generates the vertices for JavaView: 369## The coloring is calculates separately. 370 ############ 371 372 javaviewJustPointBlock:=function(vertices,group,precision) 373 local vector2floatString, cd, coordstrings, pstring; 374 375 vector2floatString:=function(v,precision) 376 local coordstring; 377 coordstring:= List(v,q->fraction2floatString(q,precision)); 378 return JoinStringsWithSeparator(coordstring," "); 379 end; 380 381 cd:=CholeskyDecomposition(GramianOfAverageScalarProductFromFiniteMatrixGroup(PointGroup(group))); 382 coordstrings:=List(vertices,i->vector2floatString(i*cd,precision)); 383 pstring:=Concatenation([List(coordstrings,s->enclosedByTag(s,"p")), 384 [enclosedByTag("3","thickness")], 385 ["\n"]]); 386 pstring:=Concatenation(pstring); 387 return enclosedByTag(pstring,"points"); 388end; 389 390 391 javaviewPointColors:=function(vertices,group) 392 local porbits, pcolors, pcolorlist; 393 394 porbits:=vertexOrbitDecomposition(vertices,group); 395 396 pcolors:=ncolorStrings(Size(porbits)); 397 pcolorlist:=List([1..Size(vertices)],v->pcolors[PositionProperty(porbits,o->v in o)]); 398 return enclosedByTag( 399 Concatenation(List(pcolorlist,s->enclosedByTag(s,"c"))), 400 "colors"); 401 402end; 403 404 405############################################################################# 406## And the edges 407################ 408 409 javaviewEdgeBlock:=function(poly,group) 410 local int2vec, vec2int, vertices, facelattice, edges, 411 edgeorbits, orientededges, orbit, firstedge, edge, 412 vecedge, map, edgeimage, estring, ecolors, 413 ecolorlist, ecstring; 414 415 int2vec:=function(n) 416 return vertices[n]; 417 end; 418 vec2int:=function(v) 419 return Position(vertices,v); 420 end; 421 vertices:=Polymake(poly,"VERTICES"); 422 facelattice:=Polymake(poly,"FACE_LATTICE"); 423 edges:=Set(facelattice[Size(facelattice)-1]); 424 edgeorbits:=edgeOrbitDecomposition(edges,vertices,group); 425 426 orientededges:=[]; 427 for orbit in edgeorbits 428 do 429 firstedge:=Set(orbit[1],int2vec); 430 Add(orientededges,List(firstedge,vec2int)); 431 for edge in orbit{[2..Size(orbit)]} 432 do 433 vecedge:=Set(edge,int2vec); 434 map:=RepresentativeActionOnRightOnSets(group,firstedge,vecedge); 435 edgeimage:=List(firstedge,v->(Concatenation(v,[1])*map)); 436 Add(orientededges,List(edgeimage,i->vec2int(i{[1..3]}))); 437 od; 438 od; 439 440 orientededges:=List(edges,e->First(orientededges,i->Set(i)=Set(e))); 441 442 estring:=List(orientededges,e->List(e-1,String)); 443 Apply(estring,s->JoinStringsWithSeparator(s," ")); 444 Apply(estring,s->enclosedByTag(s,"l")); 445 estring:=enclosedByTag( 446 Concatenation(Concatenation(estring), 447 enclosedByTag("4","thickness")) 448 ,"lines"); 449 450 ecolors:=ncolorStrings(Size(edgeorbits)); 451 ecolorlist:=List(edges,v->ecolors[PositionProperty(edgeorbits,o->v in o)]); 452 ecstring:=enclosedByTag( 453 Concatenation(List(ecolorlist,s->enclosedByTag(s,"c"))), 454 "colors"); 455 456 return enclosedByTagWithOptions(Concatenation(["\n",estring,ecstring]), 457 "lineSet", "line=\"show\" color=\"show\" arrow=\"show\""); 458 end; 459 460############################################################################# 461## FACETS 462 463 464javaviewFacetBlock:=function(poly,group) 465 local crossProduct, facets, facelattice, edges, ofacets, 466 facet, facetedges, ofacet, nextvertex, nextedge, 467 vertices, pos, v1, v2, normal, eq, testpoint, 468 fstring, facetorbits, fcolors, fcolorlist, fcstring; 469 470 crossProduct:=function(x,y) 471 return [ x[2]*y[3]-x[3]*y[2], 472 -(x[1]*y[3]-x[3]*y[1]), 473 x[1]*y[2]-x[2]*y[1]]; 474 end; 475 476 Polymake(poly,"VERTICES_IN_FACETS FACE_LATTICE"); 477 facets:=Polymake(poly,"VERTICES_IN_FACETS"); 478 ## first, generate "ordered" versions of the facets: 479 facelattice:=Polymake(poly,"FACE_LATTICE"); 480 edges:=Set(facelattice[Size(facelattice)-1]); 481 ofacets:=[]; 482 for facet in facets 483 do 484 facetedges:=Set(Filtered(edges,e->IsSubset(facet,e))); 485 ofacet:=List(facetedges[1]); 486 RemoveSet(facetedges,ofacet); 487 nextvertex:=ofacet[2]; 488 while facetedges<>[] 489 do 490 nextedge:=First(facetedges,f->nextvertex in f); 491 RemoveSet(facetedges,nextedge); 492 nextvertex:=First(nextedge,i->i<>nextvertex); 493 Add(ofacet,nextvertex); 494 od; 495 Remove(ofacet); 496 Add(ofacets,ofacet); 497 od; 498 499 ## now check if the ordering must be reversed: 500 ## We're lucky this happens in 3-space. 501 502 vertices:=Polymake(poly,"VERTICES"); 503 for pos in [1..Size(ofacets)] 504 do 505 v1:=vertices[ofacets[pos][2]]-vertices[ofacets[pos][1]]; 506 v2:=vertices[ofacets[pos][3]]-vertices[ofacets[pos][1]]; 507 normal:=crossProduct(v1,v2); 508 eq:=Concatenation([vertices[ofacets[pos][1]]*normal],normal); 509 testpoint:=vertices[Representative(Difference([1..Size(vertices)],ofacets[pos]))]; 510 if WhichSideOfHyperplane(testpoint,eq)=-1 511 then 512 ofacets[pos]:=Reversed(ofacets[pos]); 513 fi; 514 od; 515 516 fstring:=List(ofacets,e->List(e-1,String)); 517 Apply(fstring,s->JoinStringsWithSeparator(s," ")); 518 Apply(fstring,s->enclosedByTag(s,"f")); 519 fstring:=enclosedByTag(Concatenation(fstring) ,"faces"); 520 facetorbits:=edgeOrbitDecomposition(facets,vertices,group); 521 522 fcolors:=ncolorStringsPale(Size(facetorbits)); 523 fcolorlist:=List(facets,v->fcolors[PositionProperty(facetorbits,o->v in o)]); 524 fcstring:=enclosedByTag( 525 Concatenation(List(fcolorlist,s->enclosedByTag(s,"c"))), 526 "colors"); 527 return enclosedByTagWithOptions(Concatenation(["\n",fstring,fcstring]), 528 "faceSet", "face=\"show\" edge=\"show\" color=\"show\""); 529end; 530 531 532 533 534############################################################################# 535## Finally, 536## generate a nice javaview file. 537 ## It takes a number of affine matrices and generates the 538 ## image under each of them. 539 ########################### 540 541 542javaviewDatastring:=function(poly,maps,group,precision) 543 local vertices, pstring, pcstring, edgeblock, facetblock, 544 allgeometries, i, gshow, numberstring, m, pointblock, 545 facetgeometry, edgegeometry; 546 547 vertices:=Polymake(poly,"VERTICES"); 548 pstring:=javaviewJustPointBlock(vertices,group,precision); 549 pcstring:=javaviewPointColors(vertices,group); 550 edgeblock:=javaviewEdgeBlock(poly,group); 551 facetblock:=javaviewFacetBlock(poly,group); 552 553 allgeometries:=[]; 554 555 for i in [1..Size(maps)] 556 do 557 if i=1 558 then 559 gshow:="\"show\""; 560 if Size(maps)=1 561 then 562 numberstring:=""; 563 else 564 numberstring:=" FD"; 565 fi; 566 else 567 gshow:="\"hide\""; 568 numberstring:=String(i-2); 569 fi; 570 m:=maps[i]; 571 pstring:=javaviewJustPointBlock(List(vertices,v-> 572 List(Concatenation(v,[1])*m){[1..3]}), 573 group, 574 precision); 575 pointblock:=enclosedByTagWithOptions( 576 Concatenation(["\n",pstring,pcstring]), 577 "pointSet", 578 "dim=\"3\" point=\"show\" color=\"show\" " 579 ); 580 facetgeometry:=enclosedByTagWithOptions( 581 Concatenation(["\n",pointblock,"\n",facetblock,"\n"]), 582 "geometry", 583 Concatenation("name=\"points and facets ",numberstring,"\""," visible=",gshow) 584 ); 585 edgegeometry:=enclosedByTagWithOptions( 586 Concatenation(["\n",pointblock,"\n",edgeblock,"\n"]), 587 "geometry", 588 Concatenation("name=\"points and edges ",numberstring,"\""," visible=",gshow) 589 590 ); 591 592 Append(allgeometries,Concatenation([facetgeometry,"\n",edgegeometry,"\n\n"])); 593 od; 594 return allgeometries; 595end; 596 597 598 599# ########### Here we are. The orbits are colored. 600# ## Now we mark a vertex for each facet to give the whole thing some 601# ## sort of orientation. 602# ## Facet normals are calculated as well. 603# 604# normals:=[]; 605# markingpoints:=[]; 606# marklines:=[]; 607# currentpointnr:=0; 608# 609# for facetorbit in facetorbits 610# do 611# facet:=facetorbit[1]; 612# facetpos:=Position(facets,facet); 613# v1:=vertices[ofacets[facetpos][1]]; 614# v2:=vertices[ofacets[facetpos][2]]; 615# v3:=vertices[ofacets[facetpos][3]]; 616# cp:=crossProduct((v2-v1)*cd,(v3-v1)*cd); 617# normal:=cp; 618# 619# maps:=List([1..Size(facetorbit)],i-> 620# RepresentativeActionOnRightOnSets(group, 621# Set(facet,int2vec), 622# Set(facetorbit[i],int2vec) 623# ) 624# ); 625# markedvertex:=facet[1]; 626# otherpointsformark:=Flat(Filtered(edges,e->markedvertex in e 627# and IsSubset(facet,e))); 628# otherpointsformark:=Filtered(Set(otherpointsformark),p->p<>markedvertex); 629# Apply(otherpointsformark,p->1/3*(int2vec(p)-int2vec(markedvertex))+int2vec(markedvertex)); 630# 631# 632# for mapnr in [1..Size(maps)] 633# do 634# map:=maps[mapnr]; 635# #normals first: 636# normalimage:=normal*LinearPartOfAffineMatOnRight(map)^cd; 637## normalimage:=normalimage{[1..3]}*cd; 638# normalimage:=normalimage*1/sqrt(normalimage*normalimage); 639# facetpos:=Position(facets,facetorbit[mapnr]); 640# normals[facetpos]:=enclosedByTag(vector2floatString(normalimage{[1..3]},precision),"n"); 641# #now the marked corners: 642# imagepoints:=List(otherpointsformark,p->Concatenation(p,[1])*map); 643# Apply(imagepoints,p->p{[1..3]}*cd); 644# 645# for point in imagepoints 646# do 647# Add(markingpoints,enclosedByTag(vector2floatString(point,precision),"p")); 648# od; 649# Add(marklines,enclosedByTag(Concatenation([String(currentpointnr)," ",String(currentpointnr+1)]),"l")); 650# currentpointnr:=currentpointnr+2; 651# od; 652# od; 653# 654# normals:=enclosedByTag(Concatenation(normals),"normals"); 655# 656# Add(marklines,Concatenation([enclosedByTag("3","thickness"),"\n"])); 657# marklines:=enclosedByTag(Concatenation(marklines),"lines"); 658# markingpoints:=enclosedByTag(Concatenation(markingpoints),"points"); 659# Add(markingpoints,'\n'); 660# 661# markblock:=Concatenation( 662# ["\n", 663# enclosedByTagWithOptions(markingpoints, 664# "pointSet", 665# "dim=\"3\" point=\"hide\""), 666# "\n", 667# enclosedByTagWithOptions(marklines, 668# "lineSet","line=\"show\"") 669# ]); 670# 671 672 ################# 673# 674# facetgeometry:=enclosedByTagWithOptions(Concatenation(["\n",pointblock,"\n",facetblock,"\n"]),"geometry","name=\"points and facets\""); 675# edgegeometry:=enclosedByTagWithOptions(Concatenation(["\n",pointblock,"\n",edgeblock,"\n"]),"geometry","name=\"points and edges\""); 676# 677## markgeometry:=enclosedByTagWithOptions(markblock,"geometry","geometry=\"show\" name=\"marks\""); 678# 679# return Concatenation([facetgeometry,"\n",edgegeometry]);#,"\n",markgeometry]); 680 681#end; 682 683 684 685 686 687 688 689 690############################################################################# 691############################################################################# 692 693 694