File: dirichlet_to_disk.cmd

package info (click to toggle)
evolver 2.70+ds-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 17,148 kB
  • sloc: ansic: 127,395; makefile: 209; sh: 98
file content (169 lines) | stat: -rw-r--r-- 5,795 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
// dirichlet_to_disk.cmd

// Programmer: Ken Brakke, brakke@susqu.edu, http://www.susqu.edu/brakke

// For mapping simply connected regions to unit disk minimizing Dirichlet
// energy to get conformal mappings.

// Prerequisites: Dirichlet_elastic energy defined, but no others.
// Usage: Remove all constraints and boundaries.  Run to_disk or 
// to_triangle.  Then evolve.

define constraint circle_con formula: x^2 + y^2 = 1
define facet attribute form_factors real[3]

set_ff := {
  foreach facet ff do
  {  ff.form_factors[1] := ff.edge[1].length^2;
     ff.form_factors[2] := -(ff.edge[1].x*ff.edge[3].x
                        + ff.edge[1].y*ff.edge[3].y);
     ff.form_factors[3] := ff.edge[3].length^2;
  }
}

// Mapping to triangle.  Especially nice since can use the three conformal
// degrees of freedom to map particular vertices to corners, and straight
// sides should let Hessian minimize Dirichlet energy in one step.

define constraint side_1_con formula: sqrt(3)*x + y = sqrt(3);
define constraint side_2_con formula: -sqrt(3)*x + y = sqrt(3);
define constraint side_3_con formula: y = 0;

to_triangle := {
   local thise,starte,counter,sidenumber,nexte,rimcount,sidenum,lambda;
   local sidecounts,cornerx,cornery,cornerv;

   define cornerx real[4];
   cornerx := { 1.0, 0.0, -1.0, 1.0 };
   define cornery real[4];
   cornery := { 0.0, sqrt(3), 0.0, 0.0 };
   define cornerv integer[3];
   cornerv := 0;

   define sidecounts integer[3];  // how many sides mapped to each triangle side

   // Check simple connectivity
   if vertex_count - edge_count + facet_count != 1 then
   { errprintf "to_disk error: Surface not simply connected. Aborting.\n";
     return;
   };
   set_ff;  // set form factors
   set facet tension 0;

   // Go around outside edges, mapping to unit circle
   unfix vertices;
   if ( valid_element(vertex[cornerv[1]]) and
           (sum(vertex[cornerv[1]].edge,valence==1) == 2) and
        valid_element(vertex[cornerv[2]]) and
           (sum(vertex[cornerv[2]].edge,valence==1) == 2) and
        valid_element(vertex[cornerv[3]]) and
           (sum(vertex[cornerv[3]].edge,valence==1) == 2) ) then
   { foreach vertex[cornerv[1]].edge where valence == 1 do
     { thise := oid;
       break;
     };
     starte := thise;
     counter := 0;
     sidenumber := 1;
     do 
     { foreach edge[thise].vertex[2].edge eee where valence == 1 and
       eee.id != thise do
       { nexte := eee.oid; break; };
       thise := nexte;
       counter += 1;
       if ( (edge[thise].vertex[1].id == cornerv[1]) or
          (edge[thise].vertex[1].id == cornerv[2]) or
          (edge[thise].vertex[1].id == cornerv[3]) ) then
       { sidecounts[sidenumber] := counter;
         counter := 0;
         sidenumber += 1;
       };
     } while thise != starte;
   }
   else
   { printf "Legal corner vertices not given, so doing my own.\n";
     rimcount := sum(edge,valence==1);
     sidecounts[1] := rimcount idiv 3;
     sidecounts[2] := rimcount idiv 3;
     sidecounts[3] := rimcount - sidecounts[1] - sidecounts[2];
     foreach edge ee where valence ==  1 do { starte := ee.id; break; };
   };
   // Now go around moving vertices to triangle sides
   thise := starte;
   for ( sidenum := 1 ; sidenum <= 3;  sidenum += 1 )
   {
     for ( counter := 1; counter <= sidecounts[sidenum] ; counter += 1 )
     { lambda := counter/sidecounts[sidenum];
       edge[thise].vertex[1].x := (1-lambda)*cornerx[sidenum]
                                   + lambda*cornerx[sidenum+1];
       edge[thise].vertex[1].y := (1-lambda)*cornery[sidenum]
                                   + lambda*cornery[sidenum+1];
       if ( space_dimension > 2) then
         edge[thise].vertex[1].x[3] := 0; 
       if  sidenum == 1 then
         set edge[thise].vertex[1] constraint side_1_con;
       if  sidenum == 2 then
         set edge[thise].vertex[1] constraint side_2_con;
       if  sidenum == 3 then
         set edge[thise].vertex[1] constraint side_3_con;

       foreach edge[thise].vertex[2].edge eee where valence == 1 and
         eee.id != thise do
       { nexte := eee.oid; break; };
       thise := nexte;
     };
  }; 
  // Big move with rim fixed, since movable rim can have negative
  // eigenvalues.
  fix vertex where on_constraint circle_con;
  hessian;
  unfix vertex;
}
 

to_disk := {
   local thise,starte,counter,nexte,rimcount;

   // Check simple connectivity
   if vertex_count - edge_count + facet_count != 1 then
   { errprintf "to_disk error: Surface not simply connected. Aborting.\n";
     return;
   };
   set_ff;  // set form factors
   set facet tension 0;

   // Go around outside edges, mapping to unit circle
   unfix vertices;
   rimcount := sum(edge,valence==1);
   foreach edge ee where valence ==  1 do { starte := ee.id; break; };
   thise := starte;
   counter := 0;
   do 
   { 
     edge[thise].vertex[1].x := cos(counter/rimcount*2*pi); 
     edge[thise].vertex[1].y := sin(counter/rimcount*2*pi); 
     if ( space_dimension > 2) then
       edge[thise].vertex[1].x[3] := 0; 
     set edge[thise].vertex[1] constraint circle_con;
     foreach edge[thise].vertex[2].edge eee where valence == 1 and
       eee.id != thise do
     { nexte := eee.oid; break; };
     thise := nexte;
     counter += 1;
   } while thise != starte;
   
  // Big move with rim fixed, since movable rim can have negative
  // eigenvalues.
  fix vertex where on_constraint circle_con;
  hessian;
  unfix vertex;
}
 
   
   
// Prerequisites: Dirichlet energy defined, but no others.
// Usage: Remove all constraints and boundaries.  Run to_disk or to_triangl.  
// Then evolve.