nx_uv_harmonic_candidate_t280.nx source
↩ module page · 83 lines · 4556 B
1// nx_uv_harmonic_candidate_t280.nx -- Sets up a harmonic boundary condition using a rational unit circle parameterization for UV grid solving.
2import "nx_uv_gram_solver_candidate_t280.nx"
3import "nx_uv_disk_validation_candidate_t280.nx"
4// q is the successful uds_prepare result for m. Its vertex indices are
5// exactly m's cut classes; it is not a replacement source asset.
6func ulh_begin(m:*i64,q:*i64,c:i64,tolerance:f64)->*UvGramSolve{
7 if q==(0 as *i64){return 0 as *UvGramSolve}
8 if q[UV_H_NT]!=m[UV_H_NT]{return 0 as *UvGramSolve};if q[UV_H_NV]!=m[UV_H_NUV]{return 0 as *UvGramSolve}
9 let s:*UvGramSolve=uls_begin(m,c,tolerance);if s==(0 as *UvGramSolve){return s}
10 if s.result!=UV_E_NOTRUN{if s.result!=UV_OK{return s}}
11 s.operatorKind=1;s.result=UV_E_NOTRUN;s.iterations=0;s.matvecs=0;s.restarts=0
12 s.boundary=sys_mmap_try(s.n*UV_I64) as *i64
13 if (s.boundary as i64)<=0{s.boundary=0 as *i64;s.result=UV_E_ARGS;return s}
14 let index:*i64=sys_mmap_try(m[UV_H_NUV]*UV_I64) as *i64
15 let walk:*i64=sys_mmap_try(s.n*2*UV_I64) as *i64
16 if (index as i64)<=0{if (walk as i64)>0{sys_munmap_direct(walk as *u8,s.n*2*UV_I64)};s.result=UV_E_ARGS;return s}
17 if (walk as i64)<=0{sys_munmap_direct(index as *u8,m[UV_H_NUV]*UV_I64);s.result=UV_E_ARGS;return s}
18 var i:i64=0;while i<s.n{index[s.classes[i]]=i;i=i+1}
19 i=0;while i<s.n*2{s.x[i]=0.0;s.r[i]=0.0;s.z[i]=0.0;s.p[i]=0.0;s.ap[i]=0.0;s.diag[i]=0.0;i=i+1}
20 let begin:i64=m[m[UV_O_CSTART]+c];let end:i64=m[m[UV_O_CSTART]+c+1];var p:i64=begin;var bad:i64=0
21 while p<end{
22 let t:i64=m[m[UV_O_CORDER]+p];var k:i64=0
23 while k<3{
24 let a:i64=uv_ti(q,t,k);let b:i64=uv_ti(q,t,(k+1)%3);let la:i64=index[a];let lb:i64=index[b]
25 s.diag[la]=s.diag[la]+1.0;s.diag[lb]=s.diag[lb]+1.0;s.diag[s.n+la]=s.diag[s.n+la]+1.0;s.diag[s.n+lb]=s.diag[s.n+lb]+1.0
26 let edge:i64=uvi_edge_find(q,a,b)
27 if edge<0{bad=1;break}
28 let eo:i64=q[UV_O_ETAB]+edge*UV_EREC
29 if q[eo+UV_ED_T1]==0{
30 if walk[la]!=0{bad=1;break};walk[la]=lb+1;s.boundary[la]=1
31 };k=k+1
32 };if bad==1{break};p=p+1
33 }
34 sys_munmap_direct(index as *u8,m[UV_H_NUV]*UV_I64)
35 var count:i64=0;var first:i64=0-1;i=0
36 while i<s.n{if s.boundary[i]==1{count=count+1;if first<0{first=i}};i=i+1}
37 if count<3{bad=1}
38 if bad==0{
39 var current:i64=first;i=0
40 while i<count{
41 if current<0{bad=1;break};if current>=s.n{bad=1;break};if s.boundary[current]!=1{bad=1;break}
42 walk[s.n+i]=current;s.boundary[current]=2;current=walk[current]-1;i=i+1
43 }
44 if current!=first{bad=1}
45 }
46 if bad==0{
47 // Rational unit circle on four monotone quadrants. No approximate trig
48 // clamp or arbitrary geometry scale is involved; UV_PIN_SPAN is UV work units.
49 i=0;while i<count{
50 let quarter:i64=(i*4)/count;let rem:i64=(i*4)%count
51 let t:f64=__f64_from_i64(rem)/__f64_from_i64(count);let den:f64=1.0+t*t
52 let a:f64=(1.0-t*t)/den;let b:f64=2.0*t/den
53 var x:f64=a;var y:f64=b
54 if quarter==1{x=0.0-b;y=a};if quarter==2{x=0.0-a;y=0.0-b};if quarter==3{x=b;y=0.0-a}
55 let v:i64=walk[s.n+i];s.x[v]=x*__f64_from_i64(UV_PIN_SPAN);s.x[s.n+v]=y*__f64_from_i64(UV_PIN_SPAN);i=i+1
56 }
57 i=0;while i<count{
58 let a:i64=walk[s.n+i];let b:i64=walk[s.n+(i+1)%count];let cc:i64=walk[s.n+(i+2)%count]
59 let cross:f64=(s.x[b]-s.x[a])*(s.x[s.n+cc]-s.x[s.n+b])-(s.x[s.n+b]-s.x[s.n+a])*(s.x[cc]-s.x[b])
60 if uls_finite(cross)==0{bad=1;break};if cross<=0.0{bad=1;break};i=i+1
61 }
62 }
63 sys_munmap_direct(walk as *u8,s.n*2*UV_I64)
64 if bad==1{s.result=UV_E_TOPOLOGY;return s}
65 s.rhs=sys_mmap_try(s.n*2*UV_I64) as *f64
66 if (s.rhs as i64)<=0{s.rhs=0 as *f64;s.result=UV_E_ARGS;return s}
67 // Eliminate fixed boundary values once; iterative residuals compare the
68 // interior operator with this retained RHS, avoiding large cancelling terms.
69 var face:i64=0
70 while face<s.nt{var edge:i64=0;while edge<3{
71 let a:i64=s.corners[face*3+edge];let b:i64=s.corners[face*3+(edge+1)%3]
72 if s.boundary[a]==0{if s.boundary[b]>0{s.rhs[a]=s.rhs[a]+s.x[b];s.rhs[s.n+a]=s.rhs[s.n+a]+s.x[s.n+b]}}
73 if s.boundary[b]==0{if s.boundary[a]>0{s.rhs[b]=s.rhs[b]+s.x[a];s.rhs[s.n+b]=s.rhs[s.n+b]+s.x[s.n+a]}}
74 edge=edge+1};face=face+1
75 }
76 if uls_residual_apply(s)!=UV_OK{s.result=ULS_NUMERIC;return s}
77 s.rhsNorm2=uls_dot(s,s.rhs,s.rhs);s.initialNorm2=uls_dot(s,s.ap,s.ap);s.trueNorm2=s.initialNorm2
78 i=0;while i<s.n*2{if s.diag[i]<=0.0{s.result=ULS_NUMERIC;return s};s.r[i]=0.0-s.ap[i];s.z[i]=s.r[i]/s.diag[i];s.p[i]=s.z[i];i=i+1}
79 s.rho=uls_dot(s,s.r,s.z)
80 if uls_finite(s.rhsNorm2)==0{s.result=ULS_NUMERIC;return s}
81 if s.rhsNorm2==0.0{s.result=UV_OK;return s}
82 if s.rho<=0.0{s.result=ULS_NUMERIC};return s
83}