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 //----------------------------------------------------------------------------
MinIndex(const Handle (FEmTool_HAssemblyTable)& Table)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 //----------------------------------------------------------------------------
MaxIndex(const Handle (FEmTool_HAssemblyTable)& Table)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 //=======================================================================
FEmTool_Assembly(const TColStd_Array2OfInteger & Dependence,const Handle (FEmTool_HAssemblyTable)& Table)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
NullifyMatrix()127 void FEmTool_Assembly::NullifyMatrix()
128 {
129 H->Init(0.);
130 IsSolved = Standard_False;
131 }
132
133
134 //=======================================================================
135 //function : AddMatrix
136 //purpose :
137 //=======================================================================
AddMatrix(const Standard_Integer Element,const Standard_Integer Dimension1,const Standard_Integer Dimension2,const math_Matrix & Mat)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 throw Standard_DomainError("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 //=======================================================================
NullifyVector()174 void FEmTool_Assembly::NullifyVector()
175 {
176 B.Init(0.);
177 }
178
179 //=======================================================================
180 //function : AddVector
181 //purpose :
182 //=======================================================================
AddVector(const Standard_Integer Element,const Standard_Integer Dimension,const math_Vector & Vec)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 //=======================================================================
Solve()204 Standard_Boolean FEmTool_Assembly::Solve()
205 {
206 IsSolved = H->Decompose();
207
208 #ifdef OCCT_DEBUG
209 if (!IsSolved) {
210 std::cout << "Solve Echec H = " << std::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 // std::cout << "FirstIndexes(" << i << ") = " << FirstIndexes(i) << std::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 //=======================================================================
Solution(math_Vector & Solution) const363 void FEmTool_Assembly::Solution(math_Vector& Solution) const
364 {
365 if(!IsSolved) throw StdFail_NotDone("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
NbGlobVar() const417 Standard_Integer FEmTool_Assembly::NbGlobVar() const
418 {
419
420 return B.Length();
421
422 }
423
GetAssemblyTable(Handle (FEmTool_HAssemblyTable)& AssTable) const424 void FEmTool_Assembly::GetAssemblyTable(Handle(FEmTool_HAssemblyTable)& AssTable) const
425 {
426 AssTable = myRefTable;
427 }
428
429
ResetConstraint()430 void FEmTool_Assembly::ResetConstraint()
431 {
432 G.Clear();
433 C.Clear();
434 }
435
NullifyConstraint()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 //=======================================================================
AddConstraint(const Standard_Integer IndexofConstraint,const Standard_Integer Element,const Standard_Integer Dimension,const math_Vector & LinearForm,const Standard_Real Value)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
592