[occt.git] / src / FEmTool / FEmTool_Assembly.cxx

// Created on: 1998-11-17
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and / or modify it
// under the terms of the GNU Lesser General Public version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <FEmTool_Assembly.ixx>
#include <FEmTool_ListIteratorOfListOfVectors.hxx>
#include <FEmTool_ListOfVectors.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <math_Matrix.hxx>

//----------------------------------------------------------------------------
// Purpose - to find min index of global variables and define
//----------------------------------------------------------------------------
static Standard_Integer MinIndex(const Handle(FEmTool_HAssemblyTable)& Table)
{
  Standard_Integer dim, el, nvar, Imin;
  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);

  nvarl = T->Lower();

  Imin = T->Value(nvarl);

  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {
      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	Imin = Min(Imin, T->Value(nvar));
      }
    }
  return Imin;
}

//----------------------------------------------------------------------------
// Purpose - to find max index of global variables 
//----------------------------------------------------------------------------
static Standard_Integer MaxIndex(const Handle(FEmTool_HAssemblyTable)& Table)
{
  Standard_Integer dim, el, nvar, Imax;
  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T = Table->Value(diml, ell);

  nvarl = T->Lower();

  Imax = T->Value(nvarl);

  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {
      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	Imax = Max(Imax, T->Value(nvar));
      }
    }
  return Imax;
}


//=======================================================================
//function : FEmTool_Assembly
//purpose  : 
//=======================================================================
FEmTool_Assembly::FEmTool_Assembly(const TColStd_Array2OfInteger& Dependence,
				   const Handle(FEmTool_HAssemblyTable)& Table):
     myDepTable(1, Dependence.ColLength(), 1, Dependence.RowLength()),
     B(MinIndex(Table), MaxIndex(Table))
     
{
  IsSolved = Standard_False;
  myDepTable = Dependence;
  myRefTable = Table;

  TColStd_Array1OfInteger FirstIndexes(1, B.Length()); FirstIndexes.Init(B.Length());

  Standard_Integer dim, el, nvar, Imin, I0 = 1 - B.Lower(), i;

  Standard_Integer diml = Table->LowerRow(), dimu = Table->UpperRow(),
                   ell = Table->LowerCol(), elu = Table->UpperCol(), nvarl, nvaru;

  Handle(TColStd_HArray1OfInteger) T;
  for(dim = diml; dim <= dimu; dim++)
    for(el = ell; el <= elu; el++) {

      T = Table->Value(dim, el);
      nvarl = T->Lower(); nvaru = T->Upper();
      Imin = T->Value(nvarl) + I0;
     
      for(nvar = nvarl; nvar <= nvaru; nvar++) 
	Imin = Min(Imin, T->Value(nvar) + I0);
      
      for(nvar = nvarl; nvar <= nvaru; nvar++) {
	i = T->Value(nvar) + I0;
	FirstIndexes(i) = Min(FirstIndexes(i), Imin);
      }
    }


  H = new FEmTool_ProfileMatrix(FirstIndexes);

  NullifyMatrix();
  NullifyVector();
}

void FEmTool_Assembly::NullifyMatrix() 
{
  H->Init(0.);
  IsSolved = Standard_False;
}


//=======================================================================
//function : AddMatrix
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddMatrix(const Standard_Integer Element,
				 const Standard_Integer Dimension1,
				 const Standard_Integer Dimension2,
				 const math_Matrix& Mat)
{

  if(myDepTable(Dimension1, Dimension2) == 0) 
    Standard_DomainError::Raise("FEmTool_Assembly::AddMatrix");

  const TColStd_Array1OfInteger & T1 = myRefTable->Value(Dimension1,Element)->Array1();
  const TColStd_Array1OfInteger & T2 = myRefTable->Value(Dimension2,Element)->Array1();

  Standard_Integer nvarl = T1.Lower(), nvaru = Min(T1.Upper(), nvarl + Mat.RowNumber() - 1);


  Standard_Integer I, J, I0 = 1 - B.Lower(), i, ii, j,

                   i0 = Mat.LowerRow() - nvarl, j0 = Mat.LowerCol() - nvarl;

  for(i = nvarl; i <= nvaru; i++) {
    I = T1(i) + I0;
    ii = i0+i;
    for(j = nvarl; j <= i; j++) {
      J = T2(j) + I0;
      H->ChangeValue(I, J) += Mat(ii, j0+j);
    }
  }

  IsSolved = Standard_False;
}


//=======================================================================
//function : NullifyVector
//purpose  : 
//=======================================================================
void FEmTool_Assembly::NullifyVector() 
{
  B.Init(0.);
}

//=======================================================================
//function : AddVector
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddVector(const Standard_Integer Element,
				 const Standard_Integer Dimension,
				 const math_Vector& Vec) 
{
  const TColStd_Array1OfInteger & T = myRefTable->Value(Dimension,Element)->Array1();
  Standard_Integer nvarl = T.Lower(), nvaru = Min(T.Upper(), nvarl + Vec.Length() - 1),
                   i0 = Vec.Lower() - nvarl;

//  Standard_Integer I, i;
  Standard_Integer i;
  
  for(i = nvarl; i <= nvaru; i++) 
    B(T(i)) += Vec(i0 + i);
  
}


//=======================================================================
//function : Solve
//purpose  : 
//=======================================================================
Standard_Boolean FEmTool_Assembly::Solve() 
{
  IsSolved = H->Decompose();

#if DEB
  if (!IsSolved) {
    cout << "Solve Echec  H = " << endl;
    H->OutM();
  }
#endif

/* ????
  while(!IsSolved && count < 5) {
    Standard_Integer i;
    for(i = 1; i <= H->RowNumber(); i++) {
      H->ChangeValue(i,i) *= 2.;
    }
    IsSolved = H->Decompose();
    count++;
  }
*/    

// Standard_Integer count = 0;
  if(!G.IsEmpty() && IsSolved) {
    // calculating H-1 Gt
    math_Vector gi(B.Lower(), B.Upper()), qi(B.Lower(), B.Upper());

    if(GHGt.IsNull() || GHGt->RowNumber() != G.Length()) {
      TColStd_Array1OfInteger FirstIndexes(1, G.Length());
//-----------------------------------------------------------------
      TColStd_Array2OfInteger H1(1, NbGlobVar(), 1, NbGlobVar()); H1.Init(1);
      Standard_Integer i, j, k, l, BlockBeg = 1, BlockEnd;
      Standard_Boolean Block, Zero;
      for(i = 2; i <= NbGlobVar(); i++) {
	BlockEnd = i - 1;
	if(!H->IsInProfile(i, BlockEnd)) {
	  // Maybe, begin of block
	  Block = Standard_True;
	  for(j = i + 1; j <= NbGlobVar(); j++) {
	    if(H->IsInProfile(j, BlockEnd)) {
	      Block = Standard_False;
	      break;
	    }
	  }
	  if(Block) {
	    for(j = i; j <= NbGlobVar(); j++) {
	      for(k = BlockBeg; k <= BlockEnd; k++) {
		H1(j, k) = 0; H1(k, j) = 0;
	      } 
	    }
	    BlockBeg = BlockEnd + 1;
	  }
	  else i = j;
	}
      }
      
      FEmTool_ListIteratorOfListOfVectors Iter1;
      FEmTool_ListIteratorOfListOfVectors Iter2;
      for(i = 1; i <= G.Length(); i++) {
	const FEmTool_ListOfVectors& Gi = G.Value(i);
	for(j = 1; j <= i; j++) {
	  const FEmTool_ListOfVectors& Gj = G.Value(j);
	  Zero = Standard_True;
	  for(Iter1.Initialize(Gi); Iter1.More(); Iter1.Next()) {
	    const Handle(TColStd_HArray1OfReal)& a = Iter1.Value();
	    for(k = a->Lower(); k <= a->Upper(); k++) {
	      for(Iter2.Initialize(Gj); Iter2.More(); Iter2.Next()) {
		const Handle(TColStd_HArray1OfReal)& b = Iter2.Value();
		for(l = b->Lower(); l <= b->Upper(); l++) {
		  if(H1(k, l) != 0) {
		    Zero = Standard_False;
		    break;
		  }
		}
		if(!Zero) break;
	      }
	      if(!Zero) break;
	    }
	    if(!Zero) break;
	  }
	  if(!Zero) {
	    FirstIndexes(i) = j;
	    break;
	  }
	}
      }
//-----------------------------------------------------------------------
//      for(i = FirstIndexes.Lower(); i <= FirstIndexes.Upper(); i++)
//	cout << "FirstIndexes(" << i << ") = " << FirstIndexes(i) << endl;
      //      FirstIndexes.Init(1); // temporary GHGt is full matrix
      GHGt = new FEmTool_ProfileMatrix(FirstIndexes);
    }

    GHGt->Init(0.);
    Standard_Integer i, j, k;
    FEmTool_ListIteratorOfListOfVectors Iter;
    
    for(i = 1; i <= G.Length(); i++) {

      const FEmTool_ListOfVectors& L = G.Value(i);
      gi.Init(0.);
      // preparing i-th line of G (or column of Gt) 
      for(Iter.Initialize(L); Iter.More(); Iter.Next()) {

	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();

	for(j = a->Lower(); j <= a->Upper(); j++) gi(j) = a->Value(j); // gi - full line of G 

      }
                        //                                     -1 t
      H->Solve(gi, qi); // solving H*qi = gi, qi is column of H  G


      //                               -1 t
      // Calculation of product M = G H G
      // for each i all elements of i-th column of M are calculated for k >= i
      for(k = i; k <= G.Length(); k++) {

	if(GHGt->IsInProfile(k, i)) {
	  Standard_Real m = 0.; // m = M(k,i)
	
	  const FEmTool_ListOfVectors& L = G.Value(k);
	
	  for(Iter.Initialize(L); Iter.More(); Iter.Next()) {

	    const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	    for(j = a->Lower(); j <= a->Upper(); j++) m += qi(j) * a->Value(j); // scalar product of
	                                             // k-th line of G and i-th column of H-1 Gt 

	  }
	
	  GHGt->ChangeValue(k, i) = m;
	}
      }

    }


    IsSolved = GHGt->Decompose();
/*    count = 0;
    while(!IsSolved && count < 5) {
      for(i = 1; i <= GHGt->RowNumber(); i++) {
	GHGt->ChangeValue(i,i) *= 2.;
      }
      IsSolved = GHGt->Decompose();
      count++;
    }*/

  }

  return IsSolved;

}


//=======================================================================
//function : Solution
//purpose  : 
//=======================================================================
void FEmTool_Assembly::Solution(math_Vector& Solution) const
{
  if(!IsSolved) StdFail_NotDone::Raise("FEmTool_Assembly::Solution");

  if(G.IsEmpty()) H->Solve(B, Solution);
  else {

    math_Vector v1(B.Lower(), B.Upper());
    H->Solve(B, v1);
    
    math_Vector l(1, G.Length()), v2(1, G.Length());
    Standard_Integer i, j;
    FEmTool_ListIteratorOfListOfVectors Iter;
    
    for(i = 1; i <= G.Length(); i++) {
      
      const FEmTool_ListOfVectors& L = G.Value(i);
      Standard_Real m = 0.;
      
      for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
	
	const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	for(j = a->Lower(); j <= a->Upper(); j++) 
	  m += v1(j) * a->Value(j); // scalar product
	// G v1
      }
      
      v2(i) = m - C.Value(i);
    }
    
    GHGt->Solve(v2, l); // Solving M*l = v2
    
    v1 = B;
    // Calculation v1 = B-Gt*l
    // v1(j) = B(j) - Gt(j,i)*l(i) = B(j) - G(i,j)*l(i)
    // 
      for(i = 1; i <= G.Length(); i++) {
	
	const FEmTool_ListOfVectors& L = G.Value(i);
	
	for(Iter.Initialize(L); Iter.More(); Iter.Next()) {
	  
	  const Handle(TColStd_HArray1OfReal)& a = Iter.Value();
	  for(j = a->Lower(); j <= a->Upper(); j++) v1(j) -= l(i) * a->Value(j);
	  
	}
	
      }
      
    H->Solve(v1, Solution);
  }
  
}

Standard_Integer FEmTool_Assembly::NbGlobVar() const
{

  return B.Length();

}

void FEmTool_Assembly::GetAssemblyTable(Handle(FEmTool_HAssemblyTable)& AssTable) const
{
  AssTable = myRefTable;
}


void FEmTool_Assembly::ResetConstraint() 
{
  G.Clear();
  C.Clear();
}

void FEmTool_Assembly::NullifyConstraint() 
{
  FEmTool_ListIteratorOfListOfVectors Iter;
  Standard_Integer i;

  for(i = 1; i <= G.Length(); i++) {
  
    C.SetValue(i, 0.);

    for(Iter.Initialize(G.Value(i)); Iter.More(); Iter.Next()) 
      Iter.Value()->Init(0.);

  }

}


//=======================================================================
//function : AddConstraint
//purpose  : 
//=======================================================================
void FEmTool_Assembly::AddConstraint(const Standard_Integer IndexofConstraint,
				     const Standard_Integer Element,
				     const Standard_Integer Dimension,
				     const math_Vector& LinearForm,
				     const Standard_Real Value) 
{
  while(G.Length() < IndexofConstraint) {
    // Add new lines in G
    FEmTool_ListOfVectors L;
    G.Append(L);
    C.Append(0.);
  }

  FEmTool_ListOfVectors& L = G.ChangeValue(IndexofConstraint);
  
  Handle(TColStd_HArray1OfInteger) Indexes = myRefTable->Value(Dimension,Element);
  Standard_Integer i, Imax = 0, Imin = NbGlobVar();

  for(i = Indexes->Lower(); i <= Indexes->Upper(); i++) {
    Imin = Min(Imin, Indexes->Value(i));
    Imax = Max(Imax, Indexes->Value(i));
  }

  Handle(TColStd_HArray1OfReal) Coeff;

  if(L.IsEmpty()) {
    Coeff = new TColStd_HArray1OfReal(Imin,Imax);
    Coeff->Init(0.);
    L.Append(Coeff);
  }
  else {
    FEmTool_ListIteratorOfListOfVectors Iter(L);
    Standard_Integer i;
    Standard_Real  s1 = 0, s2 = 0;
    Handle(TColStd_HArray1OfReal) Aux1, Aux2;
    for(i=1; Iter.More(); Iter.Next(), i++) {
      if(Imin >= Iter.Value()->Lower()) {
	s1 = i;
	Aux1 = Iter.Value();
	if(Imax <= Iter.Value()->Upper()) {
	  s2 = s1;
	  Coeff = Iter.Value();
	  break;
	}
      }

      if(Imax <= Iter.Value()->Upper()) {
	s2 = i;
	Aux2 = Iter.Value();
      }
    }

    if(s1 != s2) {
      if(s1 == 0) {
	if(Imax < Aux2->Lower()) {
	  // inserting before first segment
	  Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	  Coeff->Init(0.);
	  L.Prepend(Coeff);
	}
	else {
	  // merge new and first segment
	  Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
	  for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	  for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	  L.First() = Coeff;
	}
      }
      else if(s2 == 0) {
	if(Imin > Aux1->Upper()) {
	  // append new
	  Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	  Coeff->Init(0.);
	  L.Append(Coeff);
	}
	else {
	  // merge new and last segment
	  Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
	  for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	  for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
	  L.Last() = Coeff;
	}
      }
      else if(Imin <= Aux1->Upper() && Imax < Aux2->Lower()) {
	// merge s1 and new
	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Imax);
	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	for(i = Aux1->Upper() + 1; i <= Imax; i++) Coeff->SetValue(i, 0.);
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	Iter.Value() = Coeff;
      }
      else if(Imin > Aux1->Upper() && Imax >= Aux2->Lower()) { 
	// merge new and first segment
	Coeff = new TColStd_HArray1OfReal(Imin, Aux2->Upper());
	for(i = Imin; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	Iter.Initialize(L);
	for(i = 1; i < s2; Iter.Next(), i++);
	Iter.Value() = Coeff;
     }
      else if(Imin > Aux1->Upper() && Imax < Aux2->Lower()) {
	// inserting new between s1 and s2
	Coeff = new TColStd_HArray1OfReal(Imin,Imax);
	Coeff->Init(0.);
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	L.InsertAfter(Coeff,Iter);
      }
      else {
	// merge s1, new, s2 and remove s2 
	Coeff = new TColStd_HArray1OfReal(Aux1->Lower(), Aux2->Upper());
	for(i = Aux1->Lower(); i <= Aux1->Upper(); i++) Coeff->SetValue(i, Aux1->Value(i));
	for(i = Aux1->Upper() + 1; i <= Aux2->Lower() - 1; i++) Coeff->SetValue(i, 0.);
	for(i = Aux2->Lower(); i <= Aux2->Upper(); i++) Coeff->SetValue(i, Aux2->Value(i));
	Iter.Initialize(L);
	for(i = 1; i < s1; Iter.Next(), i++);
	Iter.Value() = Coeff;
	Iter.Next();
	L.Remove(Iter);
      }
    }
  }

  // adding 
  Standard_Integer j = LinearForm.Lower();
  for(i = Indexes->Lower(); i <= Indexes->Upper(); i++, j++) {
    Coeff->ChangeValue(Indexes->Value(i)) += LinearForm(j);
  }
  
  C.ChangeValue(IndexofConstraint) += Value;


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