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