0026377: Passing Handle objects as arguments to functions as non-const reference...
[occt.git] / src / FEmTool / FEmTool_Assembly.cxx
1 // Created on: 1998-11-17
2 // Created by: Igor FEOKTISTOV
3 // Copyright (c) 1998-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17
18 #include <FEmTool_Assembly.hxx>
19 #include <FEmTool_ListIteratorOfListOfVectors.hxx>
20 #include <FEmTool_ListOfVectors.hxx>
21 #include <FEmTool_ProfileMatrix.hxx>
22 #include <math_Matrix.hxx>
23 #include <Standard_DimensionError.hxx>
24 #include <Standard_DomainError.hxx>
25 #include <StdFail_NotDone.hxx>
26 #include <TColStd_HArray1OfReal.hxx>
27
28 //----------------------------------------------------------------------------
29 // Purpose - to find min index of global variables and define
30 //----------------------------------------------------------------------------
31 static Standard_Integer MinIndex(const Handle(FEmTool_HAssemblyTable)& Table)
32 {
33   Standard_Integer dim, el, nvar, Imin;
34   Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
35                    ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
36
37   Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);
38
39   nvarl = T->Lower();
40
41   Imin = T->Value(nvarl);
42
43   for(dim = diml; dim <= dimu; dim++)
44     for(el = ell; el <= elu; el++) {
45       T = Table->Value(dim, el);
46       nvarl = T->Lower(); nvaru = T->Upper();
47       for(nvar = nvarl; nvar <= nvaru; nvar++) {
48         Imin = Min(Imin, T->Value(nvar));
49       }
50     }
51   return Imin;
52 }
53
54 //----------------------------------------------------------------------------
55 // Purpose - to find max index of global variables 
56 //----------------------------------------------------------------------------
57 static Standard_Integer MaxIndex(const Handle(FEmTool_HAssemblyTable)& Table)
58 {
59   Standard_Integer dim, el, nvar, Imax;
60   Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
61                    ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
62
63   Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);
64
65   nvarl = T->Lower();
66
67   Imax = T->Value(nvarl);
68
69   for(dim = diml; dim <= dimu; dim++)
70     for(el = ell; el <= elu; el++) {
71       T = Table->Value(dim, el);
72       nvarl = T->Lower(); nvaru = T->Upper();
73       for(nvar = nvarl; nvar <= nvaru; nvar++) {
74         Imax = Max(Imax, T->Value(nvar));
75       }
76     }
77   return Imax;
78 }
79
80
81
82 //=======================================================================
83 //function : FEmTool_Assembly
84 //purpose  : 
85 //=======================================================================
86 FEmTool_Assembly::FEmTool_Assembly(const TColStd_Array2OfInteger& Dependence,
87                                    const Handle(FEmTool_HAssemblyTable)& Table):
88      myDepTable(1, Dependence.ColLength(), 1, Dependence.RowLength()),
89      B(MinIndex(Table), MaxIndex(Table))
90      
91 {
92   IsSolved = Standard_False;
93   myDepTable = Dependence;
94   myRefTable = Table;
95
96   TColStd_Array1OfInteger FirstIndexes(1, B.Length()); FirstIndexes.Init(B.Length());
97
98   Standard_Integer dim, el, nvar, Imin, I0 = 1 - B.Lower(), i;
99
100   Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
101                    ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;
102
103   Handle(TColStd_HArray1OfInteger) T;
104   for(dim = diml; dim <= dimu; dim++)
105     for(el = ell; el <= elu; el++) {
106
107       T = Table->Value(dim, el);
108       nvarl = T->Lower(); nvaru = T->Upper();
109       Imin = T->Value(nvarl) + I0;
110      
111       for(nvar = nvarl; nvar <= nvaru; nvar++) 
112         Imin = Min(Imin, T->Value(nvar) + I0);
113       
114       for(nvar = nvarl; nvar <= nvaru; nvar++) {
115         i = T->Value(nvar) + I0;
116         FirstIndexes(i) = Min(FirstIndexes(i), Imin);
117       }
118     }
119
120
121   H = new FEmTool_ProfileMatrix(FirstIndexes);
122
123   NullifyMatrix();
124   NullifyVector();
125 }
126
127 void FEmTool_Assembly::NullifyMatrix() 
128 {
129   H->Init(0.);
130   IsSolved = Standard_False;
131 }
132
133
134 //=======================================================================
135 //function : AddMatrix
136 //purpose  : 
137 //=======================================================================
138 void FEmTool_Assembly::AddMatrix(const Standard_Integer Element,
139                                  const Standard_Integer Dimension1,
140                                  const Standard_Integer Dimension2,
141                                  const math_Matrix& Mat)
142 {
143
144   if(myDepTable(Dimension1, Dimension2) == 0) 
145     Standard_DomainError::Raise("FEmTool_Assembly::AddMatrix");
146
147   const TColStd_Array1OfInteger & T1 = myRefTable->Value(Dimension1,Element)->Array1();
148   const TColStd_Array1OfInteger & T2 = myRefTable->Value(Dimension2,Element)->Array1();
149
150   Standard_Integer nvarl = T1.Lower(), nvaru = Min(T1.Upper(), nvarl + Mat.RowNumber() - 1);
151
152
153   Standard_Integer I, J, I0 = 1 - B.Lower(), i, ii, j,
154
155                    i0 = Mat.LowerRow() - nvarl, j0 = Mat.LowerCol() - nvarl;
156
157   for(i = nvarl; i <= nvaru; i++) {
158     I = T1(i) + I0;
159     ii = i0+i;
160     for(j = nvarl; j <= i; j++) {
161       J = T2(j) + I0;
162       H->ChangeValue(I, J) += Mat(ii, j0+j);
163     }
164   }
165
166   IsSolved = Standard_False;
167 }
168
169
170 //=======================================================================
171 //function : NullifyVector
172 //purpose  : 
173 //=======================================================================
174 void FEmTool_Assembly::NullifyVector() 
175 {
176   B.Init(0.);
177 }
178
179 //=======================================================================
180 //function : AddVector
181 //purpose  : 
182 //=======================================================================
183 void FEmTool_Assembly::AddVector(const Standard_Integer Element,
184                                  const Standard_Integer Dimension,
185                                  const math_Vector& Vec) 
186 {
187   const TColStd_Array1OfInteger & T = myRefTable->Value(Dimension,Element)->Array1();
188   Standard_Integer nvarl = T.Lower(), nvaru = Min(T.Upper(), nvarl + Vec.Length() - 1),
189                    i0 = Vec.Lower() - nvarl;
190
191 //  Standard_Integer I, i;
192   Standard_Integer i;
193   
194   for(i = nvarl; i <= nvaru; i++) 
195     B(T(i)) += Vec(i0 + i);
196   
197 }
198
199
200 //=======================================================================
201 //function : Solve
202 //purpose  : 
203 //=======================================================================
204 Standard_Boolean FEmTool_Assembly::Solve() 
205 {
206   IsSolved = H->Decompose();
207
208 #ifdef OCCT_DEBUG
209   if (!IsSolved) {
210     cout << "Solve Echec  H = " << endl;
211     H->OutM();
212   }
213 #endif
214
215 /* ????
216   while(!IsSolved && count < 5) {
217     Standard_Integer i;
218     for(i = 1; i <= H->RowNumber(); i++) {
219       H->ChangeValue(i,i) *= 2.;
220     }
221     IsSolved = H->Decompose();
222     count++;
223   }
224 */    
225
226 // Standard_Integer count = 0;
227   if(!G.IsEmpty() && IsSolved) {
228     // calculating H-1 Gt
229     math_Vector gi(B.Lower(), B.Upper()), qi(B.Lower(), B.Upper());
230
231     if(GHGt.IsNull() || GHGt->RowNumber() != G.Length()) {
232       TColStd_Array1OfInteger FirstIndexes(1, G.Length());
233 //-----------------------------------------------------------------
234       TColStd_Array2OfInteger H1(1, NbGlobVar(), 1, NbGlobVar()); H1.Init(1);
235       Standard_Integer i, j, k, l, BlockBeg = 1, BlockEnd;
236       Standard_Boolean Block, Zero;
237       for(i = 2; i <= NbGlobVar(); i++) {
238         BlockEnd = i - 1;
239         if(!H->IsInProfile(i, BlockEnd)) {
240           // Maybe, begin of block
241           Block = Standard_True;
242           for(j = i + 1; j <= NbGlobVar(); j++) {
243             if(H->IsInProfile(j, BlockEnd)) {
244               Block = Standard_False;
245               break;
246             }
247           }
248           if(Block) {
249             for(j = i; j <= NbGlobVar(); j++) {
250               for(k = BlockBeg; k <= BlockEnd; k++) {
251                 H1(j, k) = 0; H1(k, j) = 0;
252               } 
253             }
254             BlockBeg = BlockEnd + 1;
255           }
256           else i = j;
257         }
258       }
259       
260       FEmTool_ListIteratorOfListOfVectors Iter1;
261       FEmTool_ListIteratorOfListOfVectors Iter2;
262       for(i = 1; i <= G.Length(); i++) {
263         const FEmTool_ListOfVectors& Gi = G.Value(i);
264         for(j = 1; j <= i; j++) {
265           const FEmTool_ListOfVectors& Gj = G.Value(j);
266           Zero = Standard_True;
267           for(Iter1.Initialize(Gi); Iter1.More(); Iter1.Next()) {
268             const Handle(TColStd_HArray1OfReal)& a = Iter1.Value();
269             for(k = a->Lower(); k <= a->Upper(); k++) {
270               for(Iter2.Initialize(Gj); Iter2.More(); Iter2.Next()) {
271                 const Handle(TColStd_HArray1OfReal)& b = Iter2.Value();
272                 for(l = b->Lower(); l <= b->Upper(); l++) {
273                   if(H1(k, l) != 0) {
274                     Zero = Standard_False;
275                     break;
276                   }
277                 }
278                 if(!Zero) break;
279               }
280               if(!Zero) break;
281             }
282             if(!Zero) break;
283           }
284           if(!Zero) {
285             FirstIndexes(i) = j;
286             break;
287           }
288         }
289       }
290 //-----------------------------------------------------------------------
291 //      for(i = FirstIndexes.Lower(); i <= FirstIndexes.Upper(); i++)
292 //      cout << "FirstIndexes(" << i << ") = " << FirstIndexes(i) << endl;
293       //      FirstIndexes.Init(1); // temporary GHGt is full matrix
294       GHGt = new FEmTool_ProfileMatrix(FirstIndexes);
295     }
296
297     GHGt->Init(0.);
298     Standard_Integer i, j, k;
299     FEmTool_ListIteratorOfListOfVectors Iter;
300     
301     for(i = 1; i <= G.Length(); i++) {
302
303       const FEmTool_ListOfVectors& L = G.Value(i);
304       gi.Init(0.);
305       // preparing i-th line of G (or column of Gt) 
306       for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
307
308         const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
309
310         for(j = a->Lower(); j <= a->Upper(); j++) gi(j) = a->Value(j); // gi - full line of G 
311
312       }
313                         //                                     -1 t
314       H->Solve(gi, qi); // solving H*qi = gi, qi is column of H  G
315
316
317       //                               -1 t
318       // Calculation of product M = G H G
319       // for each i all elements of i-th column of M are calculated for k >= i
320       for(k = i; k <= G.Length(); k++) {
321
322         if(GHGt->IsInProfile(k, i)) {
323           Standard_Real m = 0.; // m = M(k,i)
324         
325           const FEmTool_ListOfVectors& aL = G.Value(k);
326         
327           for(Iter.Initialize(aL); Iter.More(); Iter.Next()) {
328
329             const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
330             for(j = a->Lower(); j <= a->Upper(); j++) m += qi(j) * a->Value(j); // scalar product of
331                                                      // k-th line of G and i-th column of H-1 Gt 
332
333           }
334         
335           GHGt->ChangeValue(k, i) = m;
336         }
337       }
338
339     }
340
341
342     IsSolved = GHGt->Decompose();
343 /*    count = 0;
344     while(!IsSolved && count < 5) {
345       for(i = 1; i <= GHGt->RowNumber(); i++) {
346         GHGt->ChangeValue(i,i) *= 2.;
347       }
348       IsSolved = GHGt->Decompose();
349       count++;
350     }*/
351
352   }
353
354   return IsSolved;
355
356 }
357
358
359 //=======================================================================
360 //function : Solution
361 //purpose  : 
362 //=======================================================================
363 void FEmTool_Assembly::Solution(math_Vector& Solution) const
364 {
365   if(!IsSolved) StdFail_NotDone::Raise("FEmTool_Assembly::Solution");
366
367   if(G.IsEmpty()) H->Solve(B, Solution);
368   else {
369
370     math_Vector v1(B.Lower(), B.Upper());
371     H->Solve(B, v1);
372     
373     math_Vector l(1, G.Length()), v2(1, G.Length());
374     Standard_Integer i, j;
375     FEmTool_ListIteratorOfListOfVectors Iter;
376     
377     for(i = 1; i <= G.Length(); i++) {
378       
379       const FEmTool_ListOfVectors& L = G.Value(i);
380       Standard_Real m = 0.;
381       
382       for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
383         
384         const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
385         for(j = a->Lower(); j <= a->Upper(); j++) 
386           m += v1(j) * a->Value(j); // scalar product
387         // G v1
388       }
389       
390       v2(i) = m - C.Value(i);
391     }
392     
393     GHGt->Solve(v2, l); // Solving M*l = v2
394     
395     v1 = B;
396     // Calculation v1 = B-Gt*l
397     // v1(j) = B(j) - Gt(j,i)*l(i) = B(j) - G(i,j)*l(i)
398     // 
399       for(i = 1; i <= G.Length(); i++) {
400         
401         const FEmTool_ListOfVectors& L = G.Value(i);
402         
403         for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
404           
405           const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
406           for(j = a->Lower(); j <= a->Upper(); j++) v1(j) -= l(i) * a->Value(j);
407           
408         }
409         
410       }
411       
412     H->Solve(v1, Solution);
413   }
414   
415 }
416
417 Standard_Integer FEmTool_Assembly::NbGlobVar() const
418 {
419
420   return B.Length();
421
422 }
423
424 void FEmTool_Assembly::GetAssemblyTable(Handle(FEmTool_HAssemblyTable)& AssTable) const
425 {
426   AssTable = myRefTable;
427 }
428
429
430 void FEmTool_Assembly::ResetConstraint() 
431 {
432   G.Clear();
433   C.Clear();
434 }
435
436 void FEmTool_Assembly::NullifyConstraint() 
437 {
438   FEmTool_ListIteratorOfListOfVectors Iter;
439   Standard_Integer i;
440
441   for(i = 1; i <= G.Length(); i++) {
442   
443     C.SetValue(i, 0.);
444
445     for(Iter.Initialize(G.Value(i)); Iter.More(); Iter.Next()) 
446       Iter.Value()->Init(0.);
447
448   }
449
450 }
451
452
453 //=======================================================================
454 //function : AddConstraint
455 //purpose  : 
456 //=======================================================================
457 void FEmTool_Assembly::AddConstraint(const Standard_Integer IndexofConstraint,
458                                      const Standard_Integer Element,
459                                      const Standard_Integer Dimension,
460                                      const math_Vector& LinearForm,
461                                      const Standard_Real Value) 
462 {
463   while(G.Length() < IndexofConstraint) {
464     // Add new lines in G
465     FEmTool_ListOfVectors L;
466     G.Append(L);
467     C.Append(0.);
468   }
469
470   FEmTool_ListOfVectors& L = G.ChangeValue(IndexofConstraint);
471   
472   Handle(TColStd_HArray1OfInteger) Indexes = myRefTable->Value(Dimension,Element);
473   Standard_Integer i, Imax = 0, Imin = NbGlobVar();
474
475   for(i = Indexes->Lower(); i <= Indexes->Upper(); i++) {
476     Imin = Min(Imin, Indexes->Value(i));
477     Imax = Max(Imax, Indexes->Value(i));
478   }
479
480   Handle(TColStd_HArray1OfReal) Coeff;
481
482   if(L.IsEmpty()) {
483     Coeff = new TColStd_HArray1OfReal(Imin,Imax);
484     Coeff->Init(0.);
485     L.Append(Coeff);
486   }
487   else {
488     FEmTool_ListIteratorOfListOfVectors Iter(L);
489     Standard_Real  s1 = 0, s2 = 0;
490     Handle(TColStd_HArray1OfReal) Aux1, Aux2;
491     for(i=1; Iter.More(); Iter.Next(), i++) {
492       if(Imin >= Iter.Value()->Lower()) {
493         s1 = i;
494         Aux1 = Iter.Value();
495         if(Imax <= Iter.Value()->Upper()) {
496           s2 = s1;
497           Coeff = Iter.Value();
498           break;
499         }
500       }
501
502       if(Imax <= Iter.Value()->Upper()) {
503         s2 = i;
504         Aux2 = Iter.Value();
505       }
506     }
507
508     if(s1 != s2) {
509       if(s1 == 0) {
510         if(Imax < Aux2->Lower()) {
511           // inserting before first segment
512           Coeff = new TColStd_HArray1OfReal(Imin,Imax);
513           Coeff->Init(0.);
514           L.Prepend(Coeff);
515         }
516         else {
517           // merge new and first segment
518           Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
519           for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
520           for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
521           L.First() = Coeff;
522         }
523       }
524       else if(s2 == 0) {
525         if(Imin > Aux1->Upper()) {
526           // append new
527           Coeff = new TColStd_HArray1OfReal(Imin,Imax);
528           Coeff->Init(0.);
529           L.Append(Coeff);
530         }
531         else {
532           // merge new and last segment
533           Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
534           for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
535           for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
536           L.Last() = Coeff;
537         }
538       }
539       else if(Imin <= Aux1->Upper() && Imax < Aux2->Lower()) {
540         // merge s1 and new
541         Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
542         for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
543         for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
544         Iter.Initialize(L);
545         for(i = 1; i < s1; Iter.Next(), i++);
546         Iter.Value() = Coeff;
547       }
548       else if(Imin > Aux1->Upper() && Imax >= Aux2->Lower()) { 
549         // merge new and first segment
550         Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
551         for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
552         for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
553         Iter.Initialize(L);
554         for(i = 1; i < s2; Iter.Next(), i++);
555         Iter.Value() = Coeff;
556      }
557       else if(Imin > Aux1->Upper() && Imax < Aux2->Lower()) {
558         // inserting new between s1 and s2
559         Coeff = new TColStd_HArray1OfReal(Imin,Imax);
560         Coeff->Init(0.);
561         Iter.Initialize(L);
562         for(i = 1; i < s1; Iter.Next(), i++);
563         L.InsertAfter(Coeff,Iter);
564       }
565       else {
566         // merge s1, new, s2 and remove s2 
567         Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Aux2->Upper());
568         for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
569         for(i = Aux1->Upper() + 1; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
570         for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
571         Iter.Initialize(L);
572         for(i = 1; i < s1; Iter.Next(), i++);
573         Iter.Value() = Coeff;
574         Iter.Next();
575         L.Remove(Iter);
576       }
577     }
578   }
579
580   // adding 
581   Standard_Integer j = LinearForm.Lower();
582   for(i = Indexes->Lower(); i <= Indexes->Upper(); i++, j++) {
583     Coeff->ChangeValue(Indexes->Value(i)) += LinearForm(j);
584   }
585   
586   C.ChangeValue(IndexofConstraint) += Value;
587
588
589
590 }
591