source: tol/tol_pkg/MulVarFunAprox/DesignC2.Monomial.Spheric.Grid.tol @ 3f198cc

1941-docker-image-to-testdevel
Last change on this file since 3f198cc was 3f198cc, checked in by Víctor de Buen <vdebuen@…>, 17 years ago

Preparando para permitir radio variante en los entornos esféricos

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1885 be4846ae-f051-483a-831f-b64e7d08527f
  • Property mode set to 100644
File size: 8.8 KB
Line 
1/* -*- mode: C++ -*- */
2//////////////////////////////////////////////////////////////////////////////
3// FILE   : DesignC2.tol
4// PURPOSE: Defines Class @DesignC2
5//////////////////////////////////////////////////////////////////////////////
6
7
8//////////////////////////////////////////////////////////////////////////////
9Text _.autodoc.member.DesignC2.Monomial.Spheric.Grid =
10"";
11Class @DesignC2.Monomial.Spheric.Grid :
12  @DesignC2.Monomial,
13  @DesignC2.Spheric
14//////////////////////////////////////////////////////////////////////////////
15{
16  @IndexedGrid.Border _.grid;
17//////////////////////////////////////////////////////////////////////////////
18//Configuration parameters
19//////////////////////////////////////////////////////////////////////////////
20  Real _.radius;
21  Real _.radius.margin;
22  Real _.grid.size;
23
24//////////////////////////////////////////////////////////////////////////////
25//Auxiliar members
26//////////////////////////////////////////////////////////////////////////////
27  VMatrix _.Y.h = Rand(0,0,0,0);
28  VMatrix _.Y.neighbourhood = Rand(0,0,0,0);
29
30  ////////////////////////////////////////////////////////////////////////////
31  Static @DesignC2.Monomial.Spheric.Grid New(
32    VMatrix X,
33    Real radius.margin,
34    Real radius,
35    Real do.build.border)
36  ////////////////////////////////////////////////////////////////////////////
37  {
38    @DesignC2.Monomial.Spheric.Grid new =
39    [[
40      VMatrix _.X = X;
41      Real    _.n = VColumns(X);
42      Real    _.J = VRows(X);
43      Real    _.radius = radius;
44      Real    _.radius.margin = radius.margin;
45      Real    _.grid.size = radius/(1+radius.margin);
46      @IndexedGrid.Border _.grid = @IndexedGrid.Border::New(
47         X,_.grid.size,True,do.build.border,do.build.border)
48    ]]
49  };
50
51  //////////////////////////////////////////////////////////////////////////////
52  VMatrix get.radius(VMatrix Y.idx)
53  //////////////////////////////////////////////////////////////////////////////
54  {
55    Constant(VRows(Y.idx),1,_.radius)
56  };
57
58  //////////////////////////////////////////////////////////////////////////////
59  VMatrix closest.reference.index(VMatrix X)
60  //////////////////////////////////////////////////////////////////////////////
61  {
62    WriteLn("TRACE [closest.reference.index] 1");
63    SubRow(_.Y.neighbourhood, Col2Set(_.grid::point2index(X)))
64  };
65
66  //////////////////////////////////////////////////////////////////////////////
67  Real build.reference.points(Real void)
68  //////////////////////////////////////////////////////////////////////////////
69  {
70    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::Grid.reference.points] ";
71    VMatrix _.Y := _.grid::_.Y;
72    VMatrix _.Y.h := Mat2VMat(_.grid::_.Y.h);
73    VMatrix _.Y.neighbourhood := _.Y.h | Group("ConcatColumns",
74    For(1,2*_.n, VMatrix(Real i)
75    {
76      VMatrix aux = SubCol(_.grid::_.Y.neighbourhood.index,[[i]]);
77      IfVMat(aux,aux,_.Y.h)
78    }));
79    Set X.class = Classify(Col2Set(_.grid::point2index(_.X)),Compare);
80    VMatrix _.Yj := Mat2VMat(SetMat([[EvalSet(X.class,Card)]]),1);
81    Real _.H := _.grid::_.H
82  };
83
84  //////////////////////////////////////////////////////////////////////////////
85  Real build.stress.points(Real void)
86  //////////////////////////////////////////////////////////////////////////////
87  {
88    Text _MID = "[@DesignC2.Spheric::build.stress.points] ";
89    WriteLn("TRACE "+_MID+"1");
90    Matrix Y.neighbour.triplet = Group("ConcatRows", For(1,_.H,Matrix(Real h)
91    {
92      Set neighbours = Unique(VRow2Set(SubRow(_.Y.neighbourhood,[[h]])));
93      Set neighbours.up = Select(neighbours,Real(Real s) { s>h });
94      Matrix s = SetCol(Sort(neighbours.up,Compare));
95      Rand(Rows(s),1,h,h)|s|Rand(Rows(s),1,1,1)
96    }));
97
98    WriteLn("TRACE "+_MID+"2");
99    VMatrix Y.neighbour.1.map = Triplet(Y.neighbour.triplet,_.H,_.H);
100    VMatrix Y.neighbour.2.map = Y.neighbour.1.map*Y.neighbour.1.map;
101    VMatrix Y.neighbour.12.map = Max(Y.neighbour.1.map,Y.neighbour.2.map);
102    Matrix Y.neighbour.12.triplet = VMat2Triplet(Y.neighbour.12.map);
103    VMatrix Z.h.a = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[1]]));
104    VMatrix Z.h.b = Mat2VMat(SubCol(Y.neighbour.12.triplet,[[2]]));
105    VMatrix Y.neighbour.a = SubRow(_.Y,VCol2Set(Z.h.a));
106    VMatrix Y.neighbour.b = SubRow(_.Y,VCol2Set(Z.h.b));
107    VMatrix Z.all = (Y.neighbour.a + Y.neighbour.b)*0.5;
108
109    WriteLn("TRACE "+_MID+"3");
110    VMatrix Z.interior.a =
111    {
112      VMatrix Z.norm2 = RPow(Y.neighbour.a-Z.all,2)*Rand(_.n,1,1,1);
113      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
114      LE(Z.tao,Rand(VRows(Z.tao),1,1,1))
115    };
116    VMatrix Z.interior.b =
117    {
118      VMatrix Z.norm2 = RPow(Y.neighbour.b-Z.all,2)*Rand(_.n,1,1,1);
119      VMatrix Z.tao = RPow(Z.norm2,0.5)/_.radius;
120      LE(Z.tao,Rand(VRows(Z.tao),1,1,1))
121    };
122    VMatrix match.interior = Z.interior.a $* Z.interior.b;
123
124    Set sel.interior = MatQuery::SelectMatch(match.interior);
125    VMatrix Z.interior = SubRow(Z.all,sel.interior);
126    VMatrix Z.h.interior = SubRow(Z.h.a | Z.h.b,sel.interior);
127    Set criterium = Range(1,_.n,1);
128    Set sel.unique = Col2Set(Sort(VMat2Mat(Z.interior),criterium,True));
129    WriteLn("TRACE "+_MID+"4");
130    VMatrix _.Z := SubRow(Z.interior, sel.unique);
131    VMatrix _.Z.h := SubRow(Z.h.interior, sel.unique);
132    Real _.Tz := VRows(_.Z);
133    VMatrix _.Z.Y.neighbour :=
134    {
135      Matrix Z.t = DifEq(1/(1-B),Rand(_.Tz,1,1,1));
136      Matrix Z.th = (Z.t << Z.t) |
137       VMat2Mat((SubCol(_.Z.h,[[1]])<<SubCol(_.Z.h,[[2]])));
138      Matrix Z.Y.idx = PivotByRows(Z.th,Sort(Z.th,[[2]]));
139      VMatrix Z.Y.dif =
140       SubRow(_.Z,Col2Set(SubCol(Z.Y.idx,[[1]]))) -
141       SubRow(_.Y,Col2Set(SubCol(Z.Y.idx,[[2]])));
142        VMatrix Z.Y.norm2 = RPow(Z.Y.dif,2)*Rand(_.n,1,1,1);
143        VMatrix Z.Y.tao = RPow(Z.Y.norm2,0.5)/_.radius;
144        VMatrix Z.Y.match = LE(Z.Y.tao,Rand(VRows(Z.Y.tao),1,1,1));
145      SubRow(Mat2VMat(Z.Y.idx),MatQuery::SelectMatch(Z.Y.match))
146    };
147    VMatrix aux = MatQuery::SortAndCount(SubCol(_.Z.Y.neighbour,[[2]]));
148    VMatrix _.Yt := Triplet(VMat2Mat(SubCol(aux,[[1]])) | 
149                            Rand(VRows(aux),1,1,1) |
150                            VMat2Mat(SubCol(aux,[[2]])),_.H,1);
151    Real _.Mz := Case(
152      prior.Z.order< 0, 0,
153      prior.Z.order==0, _.Tz,
154      prior.Z.order==1, _.Tz*(1+_.n),
155      prior.Z.order==2, _.Tz*(1+_.n*3/2+(1/2)*_.n^2));
156    _.Tz
157  };
158
159  //////////////////////////////////////////////////////////////////////////////
160  Real build.border.points(Real void)
161  //////////////////////////////////////////////////////////////////////////////
162  {
163    If(!Rows(_.grid::_.border.assign),{
164      VMatrix _.Yb := Zeros(_.H,1);
165      Real _.Tb := 0
166    },{
167    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.border.points] ";
168    VMatrix   B.t  = Mat2VMat(SubCol(_.grid::_.border.assign,[[1]]));
169    VMatrix _.B.h := Mat2VMat(SubCol(_.grid::_.border.assign,[[2]]));
170    VMatrix _.B := SubRow(_.grid::_.border::_.Y, VCol2Set(B.t));
171    VMatrix aux = MatQuery::SortAndCount(_.B.h);
172    VMatrix _.Yb := Triplet(VMat2Mat(SubCol(aux,[[1]])) | 
173                            Rand(VRows(aux),1,1,1) |
174                            VMat2Mat(SubCol(aux,[[2]])),_.H,1);
175
176    Real _.Tb := VRows(_.B);
177    Real _.Mb := Case(
178      prior.B.order< 0, 0,
179      prior.B.order==0, _.Tb,
180      prior.B.order==1, _.Tb*(1+_.n),
181      prior.B.order==2, _.Tb*(1+_.n*3/2+(1/2)*_.n^2));
182    _.Tb
183    })
184  };
185
186  //////////////////////////////////////////////////////////////////////////////
187    Real build.model.design(Real void)
188  //////////////////////////////////////////////////////////////////////////////
189  {
190    Text _MID = "[@DesignC2.Monomial.Spheric.Grid::build.model.design] ";
191    FunRn2R::@BaseC2 base.max =
192      FunRn2R::@BaseC2.Monomial::New(_.n,max.deg);
193    Set Append(_.base, [[base.max]]);
194    Real _.base.sel := Card(_.base);
195    Real _.K0 := base.max::_.K;
196    _build.model.design(void)
197  };
198
199  ///////////////////////////////////////////////////////////////////////////
200    VMatrix generate.internal.points(
201      Real m,
202      Real force.reference.near.nodal)
203  ///////////////////////////////////////////////////////////////////////////
204  {
205    VMatrix h = Round(Rand(m,1,0.5,_.H+0.499999999));
206    VMatrix y = SubRow(_.Y,VCol2Set(h));
207    VMatrix v = Gaussian(m,_.n,0,1);
208    VMatrix u = v*(VMatDat(MtMSqr(v),1,1)^-0.5);
209    y+u*(Rand(0,1)*_.radius)
210  };
211
212  ////////////////////////////////////////////////////////////////////////////
213    VMatrix match.evaluable(VMatrix X)
214  ////////////////////////////////////////////////////////////////////////////
215  {
216    Real m = VRows(X);
217    Mat2VMat(_.grid::point2index(X))
218/*
219    VMatrix X.closest.index = closest.reference.index(X);
220    VMatrix X.h=SubRow(SubCol(_.Y.neighbourhood,[[1]]),
221                       VCol2Set(X.closest.index));
222    VMatrix X.dif = X - SubRow(_.Y,VCol2Set(X.h));
223    VMatrix X.norm2 = RPow(X.dif,2)*Rand(_.n,1,1,1);
224    VMatrix X.tao = RPow(X.norm2,0.5)/_.radius;
225    LE(X.tao,Rand(m,1,1,1))
226*/
227  }
228
229};
Note: See TracBrowser for help on using the repository browser.