NeoPZ
pzmatred.cpp
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <stdio.h>
8 #include <fstream>
9 using namespace std;
10 
11 
12 #include "pzmatred.h"
13 #include "pzfmatrix.h"
14 #include "pzstepsolver.h"
15 #include "tpzverysparsematrix.h"
16 
17 #include "TPZPersistenceManager.h"
18 
19 #include <sstream>
20 #include "pzlog.h"
21 #ifdef LOG4CXX
22 static LoggerPtr logger(Logger::getLogger("pz.matrix.tpzmatred"));
23 #endif
24 
25 /*************************** Public ***************************/
26 
27 /******************/
28 /*** Construtor ***/
29 
30 template<class TVar, class TSideMatrix>
31 TPZMatRed<TVar, TSideMatrix>::TPZMatRed () :
32 TPZRegisterClassId(&TPZMatRed::ClassId),
33 TPZMatrix<TVar>( 0, 0 ), fK11(0,0),fK01(0,0),fK10(0,0),fF0(0,0),fF1(0,0), fMaxRigidBodyModes(0), fNumberRigidBodyModes(0)
34 {
35  fDim0=0;
36  fDim1=0;
37  fK01IsComputed = 0;
38  fIsReduced = 0;
39  fF0IsComputed = false;
40 }
41 
42 template<class TVar, class TSideMatrix>
43 TPZMatRed<TVar, TSideMatrix>::TPZMatRed( int64_t dim, int64_t dim00 ):
44 TPZRegisterClassId(&TPZMatRed::ClassId),
45 TPZMatrix<TVar>( dim,dim ), fK11(dim-dim00,dim-dim00,0.), fK01(dim00,dim-dim00,0.),
46 fK10(dim-dim00,dim00,0.), fF0(dim00,1,0.),fF1(dim-dim00,1,0.), fMaxRigidBodyModes(0), fNumberRigidBodyModes(0)
47 {
48  if(dim<dim00) TPZMatrix<TVar>::Error(__PRETTY_FUNCTION__,"dim k00> dim");
49  fDim0=dim00;
50  fDim1=dim-dim00;
51  fK01IsComputed = 0;
52  fF0IsComputed = false;
53  fIsReduced = 0;
54 }
55 
56 template<class TVar, class TSideMatrix>
57 TPZMatRed<TVar, TSideMatrix >::~TPZMatRed(){
58 }
59 
60 template<class TVar, class TSideMatrix>
61 int TPZMatRed<TVar, TSideMatrix>::IsSimetric() const {
62  if(fK00) return this->fK00->IsSimetric();
63  return 0;
64 }
65 
66 template<class TVar, class TSideMatrix>
67 void TPZMatRed<TVar, TSideMatrix>::SimetrizeMatRed() {
68  // considering fK00 is simetric, only half of the object is assembled.
69  // this method simetrizes the matrix object
70 
71  if(!fK00 || !this->fK00->IsSimetric()) return;
72  fK01.Transpose(&fK10);
73 
74  fK11.Simetrize();
75 
76  /*
77  int64_t row,col;
78  for(row=0; row<fDim1; row++) {
79  for(col=row+1; col<fDim1; col++) {
80  (fK11)(col,row) = (fK11)(row,col);
81  }
82  }
83  */
84 }
85 
86 template<class TVar, class TSideMatrix>
87 int
88 TPZMatRed<TVar, TSideMatrix>::PutVal(const int64_t r,const int64_t c,const TVar& value ){
89  int64_t row(r),col(c);
90  if (IsSimetric() && row > col ) Swap( &row, &col );
91  if (row<fDim0 && col<fDim0) fK00->PutVal(row,col,value);
92  if (row<fDim0 && col>=fDim0) fK01.PutVal(row,col-fDim0,(TVar)value);
93  if (row>=fDim0 && col<fDim0) fK10.PutVal(row-fDim0,col,(TVar)value);
94  if (row>=fDim0 && col>=fDim0) fK11.PutVal(row-fDim0,col-fDim0,(TVar)value);
95 
96  return( 1 );
97 }
98 
99 template<class TVar, class TSideMatrix>
100 const TVar&
101 TPZMatRed<TVar,TSideMatrix>::GetVal(const int64_t r,const int64_t c ) const {
102  int64_t row(r),col(c);
103 
104  if (IsSimetric() && row > col ) Swap( &row, &col );
105  if (row<fDim0 && col<fDim0) return ( fK00->GetVal(row,col) );
106  if (row<fDim0 && col>=fDim0) return ( fK01.GetVal(row,col-fDim0) );
107  if (row>=fDim0 && col<fDim0) return ( fK10.GetVal(row-fDim0,col) );
108  return (fK11.GetVal(row-fDim0,col-fDim0) );
109 
110 }
111 
112 template<class TVar, class TSideMatrix>
113 TVar& TPZMatRed<TVar,TSideMatrix>::s(const int64_t r,const int64_t c ) {
114  int64_t row(r),col(c);
115 
116  if (r < fDim0 && IsSimetric() && row > col ) Swap( &row, &col );
117  if (row<fDim0 && col<fDim0) return ( fK00->s(row,col) );
118  if (row<fDim0 && col>=fDim0) return ( (TVar &)fK01.s(row,col-fDim0) );
119  if (row>=fDim0 && col<fDim0) return ( (TVar &)(fK10.s(row-fDim0,col)) );
120  return ((TVar &)(fK11.s(row-fDim0,col-fDim0)) );
121 
122 }
123 
124 template<class TVar, class TSideMatrix>
125 void TPZMatRed<TVar,TSideMatrix>::SetSolver(TPZAutoPointer<TPZMatrixSolver<TVar> > solver)
126 {
127  fK00=solver->Matrix();
128  fSolver = solver;
129 }
130 
131 
132 template<class TVar,class TSideMatrix>
133 void
134 TPZMatRed<TVar, TSideMatrix>::SetK00(TPZAutoPointer<TPZMatrix<TVar> > K00)
135 {
136  fK00=K00;
137 }
138 
139 template<class TVar, class TSideMatrix>
140 void TPZMatRed<TVar,TSideMatrix>::SetF(const TPZFMatrix<TVar> & F)
141 {
142 
143  int64_t FCols=F.Cols(),c,r,r1;
144 
145  fF0.Redim(fDim0,FCols);
146  fF1.Redim(fDim1,FCols);
147  fF0IsComputed = false;
148 
149  for(c=0; c<FCols; c++){
150  r1=0;
151  for(r=0; r<fDim0; r++){
152  fF0.PutVal( r,c,F.GetVal(r,c) ) ;
153  }
154  //aqui r=fDim0
155  for( ;r<fDim0+fDim1; r++){
156  fF1.PutVal( r1++,c,F.GetVal(r,c) );
157  }
158  }
159 #ifdef LOG4CXX
160  if (logger->isDebugEnabled()) {
161  std::stringstream sout;
162  F.Print("F Input",sout);
163  fF0.Print("fF0 Initialized",sout);
164  fF1.Print("fF1 Initialized",sout);
165  LOGPZ_DEBUG(logger, sout.str())
166  }
167 #endif
168 }
169 
170 template<class TVar, class TSideMatrix>
171 void TPZMatRed<TVar, TSideMatrix>::F1Red(TPZFMatrix<TVar> &F1Red)
172 {
173  if (!fDim0)
174  {
175  F1Red = fF1;
176  return;
177  }
178 #ifdef LOG4CXX
179  if (logger->isDebugEnabled()) {
180  std::stringstream sout;
181  sout << "fF0 input " << std::endl;
182  fF0.Print("fF0 = ",sout,EMathematicaInput);
183  LOGPZ_DEBUG(logger, sout.str())
184  }
185 #endif
186  F1Red.Resize(fK11.Rows(),fF0.Cols());
187  if (!fF0IsComputed)
188  {
189  DecomposeK00();
190  fSolver->Solve(fF0,fF0);
191  fF0IsComputed = true;
192  }
193 #ifdef LOG4CXX
194  if (logger->isDebugEnabled()) {
195  std::stringstream sout;
196  sout << "After computing F0Invert" << std::endl;
197  fF0.Print("F0Invert",sout,EMathematicaInput);
198  LOGPZ_DEBUG(logger, sout.str())
199  }
200 #endif
201 
202 #ifdef LOG4CXX
203  if (logger->isDebugEnabled()) {
204  std::stringstream sout;
205  sout << "Input fF1" << std::endl;
206  fF1.Print("fF1",sout);
207  LOGPZ_DEBUG(logger, sout.str())
208  }
209 #endif
210  //make [F1]=[F1]-[K10][F0Invert]
211  fK10.MultAdd((fF0),fF1,(F1Red),-1,1);
212 #ifdef LOG4CXX
213  if (logger->isDebugEnabled()) {
214  std::stringstream sout;
215  F1Red.Print("F1 Reduced", sout);
216  LOGPZ_DEBUG(logger, sout.str())
217  }
218 #endif
219  return;
220 }
221 
222 template<class TVar, class TSideMatrix>
223 void TPZMatRed<TVar,TSideMatrix>::K11Reduced(TPZFMatrix<TVar> &K11, TPZFMatrix<TVar> &F1)
224 {
225 
226  if(!fK01IsComputed)
227  {
228  DecomposeK00();
229  SimetrizeMatRed();
230  fSolver->Solve(fK01,fK01);
231  TPZStepSolver<TVar> *step = dynamic_cast<TPZStepSolver<TVar> *>(fSolver.operator->());
232  if (step->Singular().size())
233  {
234  std::cout << "Address " << (void *) step << " Number of singular modes " << step->Singular().size() << std::endl;
235  }
236 #ifdef LOG4CXX
237  if(logger->isDebugEnabled())
238  {
239  std::stringstream sout;
240  sout << "K01 after reduction\n";
241  fK01.Print("fK01 = ",sout,EMathematicaInput);
242  LOGPZ_DEBUG(logger, sout.str())
243  }
244 #endif
245  fK01IsComputed = true;
246  }
247  fK10.MultAdd(fK01,fK11,(K11),-1.,1.);
248  F1Red(F1);
249 
250  return;
251 }
252 
253 #include "tpzverysparsematrix.h"
254 
255 template<>
256 void TPZMatRed<double, TPZVerySparseMatrix<double> >::K11Reduced(TPZFMatrix<double> &K11, TPZFMatrix<double> &F1)
257 {
258  std::cout << __PRETTY_FUNCTION__ << " should never be called\n";
259  static TPZFMatrix<double> temp;
260  return;
261 }
262 
263 template<>
264 void TPZMatRed<float, TPZVerySparseMatrix<float> >::K11Reduced(TPZFMatrix<float> &K11, TPZFMatrix<float> &F1)
265 {
266  std::cout << __PRETTY_FUNCTION__ << " should never be called\n";
267  static TPZFMatrix<float> temp;
268  return;
269 }
270 
271 template<>
272 void TPZMatRed<long double, TPZVerySparseMatrix<long double> >::K11Reduced(TPZFMatrix<long double> &K11, TPZFMatrix<long double> &F1)
273 {
274  std::cout << __PRETTY_FUNCTION__ << " should never be called\n";
275  static TPZFMatrix<long double> temp;
276  return;
277 }
278 
279 template<>
280 void TPZMatRed<std::complex<double>, TPZVerySparseMatrix<std::complex<double> > >::K11Reduced(TPZFMatrix<std::complex<double> > &K11, TPZFMatrix<std::complex<double> > &F1)
281 {
282  std::cout << __PRETTY_FUNCTION__ << " should never be called\n";
283  static TPZFMatrix<std::complex<double> > temp;
284  return;
285 }
286 
287 
288 template<class TVar ,class TSideMatrix>
289 void TPZMatRed<TVar,TSideMatrix>::U1(TPZFMatrix<TVar> & F)
290 {
291  TPZFNMatrix<1000,TVar> K1Red(fDim1,fDim1), F1Red(fDim1,fF1.Cols());
292  K11Reduced(K1Red, F1Red);
293  F=(F1Red);
294  K1Red.SolveDirect( F ,ELU);
295 
296 
297 }
298 
299 template<>
300 void TPZMatRed<REAL, TPZVerySparseMatrix<REAL> >::UGlobal(const TPZFMatrix<REAL> & U1, TPZFMatrix<REAL> & result)
301 {
302  //[u0]=[A00^-1][F0]-[A00^-1][A01]
303  //compute [F0]=[A00^-1][F0]
304  if (!fF0IsComputed)
305  {
306  DecomposeK00();
307  fSolver->Solve(fF0,fF0);
308  fF0IsComputed = true;
309  }
310  if(!fK01IsComputed)
311  {
312  TPZFMatrix<REAL> k01(fK01);
313  DecomposeK00();
314  fSolver->Solve(k01,k01);
315  fK01 = k01;
316  fK01IsComputed = 1;
317  }
318 
319  //make [u0]=[F0]-[U1]
320  TPZFMatrix<REAL> u0( fF0.Rows() , fF0.Cols() );
321  fK01.MultAdd(U1,(fF0),u0,-1,0);
322 
323  result.Redim( fDim0+fDim1,fF0.Cols() );
324  int64_t c,r,r1;
325 
326  for(c=0; c<fF0.Cols(); c++)
327  {
328  r1=0;
329  for(r=0; r<fDim0; r++)
330  {
331  result.PutVal( r,c,u0.GetVal(r,c) ) ;
332  }
333  //aqui r=fDim0
334  for( ;r<fDim0+fDim1; r++)
335  {
336  result.PutVal( r,c,U1.GetVal(r1++,c) );
337  }
338  }
339 }
340 
341 template<class TVar, class TSideMatrix>
342 void TPZMatRed<TVar, TSideMatrix >::UGlobal(const TPZFMatrix<TVar> & U1, TPZFMatrix<TVar> & result)
343 {
344  TPZFMatrix<TVar> u0( fF0.Rows() , fF0.Cols() );
345 
346  if(fK01IsComputed)
347  {
348  //[u0]=[A00^-1][F0]-[A00^-1][A01]
349  //compute [F0]=[A00^-1][F0]
350  if(!fF0IsComputed)
351  {
352  DecomposeK00();
353  fSolver->Solve(fF0,fF0);
354  fF0IsComputed = true;
355  }
356  //make [u0]=[F0]-[U1]
357  fK01.MultAdd(U1,(fF0),u0,-1,1);
358  } else {
359  TPZFMatrix<TVar> K01U1(fK01.Rows(),U1.Cols(),0.);
360  fK01.MultAdd(U1,fF0,K01U1,-1.,1.);
361  DecomposeK00();
362  fSolver->Solve(K01U1, u0);
363  }
364 
365  //compute result
366 #ifdef LOG4CXX
367  if(logger->isDebugEnabled())
368  {
369  std::stringstream sout;
370  U1.Print("U1 = ",sout,EMathematicaInput);
371  fF0.Print("fF0 ",sout,EMathematicaInput);
372  u0.Print("u0 " ,sout,EMathematicaInput);
373  LOGPZ_DEBUG(logger,sout.str())
374 
375  }
376 #endif
377 
378  result.Redim( fDim0+fDim1,fF0.Cols() );
379  int64_t c,r,r1;
380 
381  for(c=0; c<fF0.Cols(); c++)
382  {
383  r1=0;
384  for(r=0; r<fDim0; r++)
385  {
386  result.PutVal( r,c,u0.GetVal(r,c) ) ;
387  }
388  //aqui r=fDim0
389  for( ;r<fDim0+fDim1; r++)
390  {
391  result.PutVal( r,c,U1.GetVal(r1++,c) );
392  }
393  }
394 }
395 
396 template<class TVar, class TSideMatrix>
397 void TPZMatRed<TVar, TSideMatrix>::UGlobal2(TPZFMatrix<TVar> & U1, TPZFMatrix<TVar> & result)
398 {
399  TPZFMatrix<TVar> u0( fF0.Rows() , fF0.Cols() );
400 
401  if(fK01IsComputed)
402  {
403  //[u0]=[A00^-1][F0]-[A00^-1][A01][u1]
404  //compute [F0]=[A00^-1][F0]
405  if(!fF0IsComputed)
406  {
407  DecomposeK00();
408  fSolver->Solve(fF0,fF0);
409  fF0IsComputed = true;
410  }
411  //make [u0]=[F0]-[U1]
412  fK01.MultAdd(U1,(fF0),u0,-1,1);
413  } else {
414  TPZFMatrix<TVar> K01U1(fK01.Rows(),U1.Cols(),0.);
415  fK01.MultAdd(U1,fF0,K01U1,-1.,1.);
416  DecomposeK00();
417  fSolver->Solve(K01U1, u0);
418  }
419 
420  //compute result
421 #ifdef LOG4CXX
422  if(logger->isDebugEnabled())
423  {
424  std::stringstream sout;
425  fF0.Print("fF0 ",sout);
426  u0.Print("u0 " ,sout);
427  LOGPZ_DEBUG(logger,sout.str())
428 
429  }
430 #endif
431 
432  result.Redim( fDim0+fDim1,fF0.Cols() );
433  int64_t nglob = fDim0+fDim1;
434 
435  int64_t c,r,r1;
436 
437  for(c=0; c<fF0.Cols(); c++)
438  {
439  r1=0;
440  for(r=0; r<fDim0; r++)
441  {
442  result(r,c) = u0(r,c) ;
443  }
444  //aqui r=fDim0
445  for( ;r<nglob; r++)
446  {
447  result(r,c) = U1(r1++,c);
448  }
449  }
450 }
451 
452 
453 template<class TVar,class TSideMatrix>
454 void TPZMatRed<TVar, TSideMatrix>::Print(const char *name , std::ostream &out ,const MatrixOutputFormat form ) const
455 {
456  if(form != EInputFormat) {
457  out << "Writing matrix 'TPZMatRed<TSideMatrix>::" << name;
458  out << "' (" << this->Dim() << " x " << this->Dim() << "):\n";
459  out << "fIsReduced " << this->fIsReduced << std::endl;
460  out << "fF0IsComputed " << this->fF0IsComputed << std::endl;
461  out << "fK01IsComputed " << this->fK01IsComputed << std::endl;
462  out << "fNumberRigidBodyModes " << this->fNumberRigidBodyModes << std::endl;
463  out << std::endl;
464  fK00->Print("K00 =",out,form);
465  fK01.Print("K01 = ",out,form);
466  fK10.Print("K10 = ",out,form);
467  fK11.Print("K11 = ",out,form);
468 
469 
470  fF0.Print("F0 = ",out,form);
471  fF1.Print("F1 = ",out,form);
472 
473  out << "Matrix norms K00 " << Norm(*fK00.operator->()) << " K01 " << Norm(fK01) << " K10 " << Norm(fK10) << " K11 " << Norm(fK11);
474  out << "\n\n";
475  } else {
476  TPZMatrix<TVar>::Print(name,out,form);
477  }
478 }
479 
480 template<class TVar, class TSideMatrix>
481 int TPZMatRed<TVar,TSideMatrix>::Redim(int64_t dim, int64_t dim00){
482  if(dim<dim00) TPZMatrix<TVar>::Error(__PRETTY_FUNCTION__,"dim k00> dim");
483  if(fK00) fK00->Redim(dim00,dim00);
484 
485  fDim0=dim00;
486  fDim1=dim-dim00;
487 
488  fK01.Redim(fDim0,fDim1);
489  fK10.Redim(fDim1,fDim0);
490  fK11.Redim(fDim1,fDim1);
491 
492  fF0.Redim(fDim0,1);
493  fF1.Redim(fDim1,1);
494  this->fRow = dim;
495  this->fCol = dim;
496  fIsReduced = false;
497  fK01IsComputed = false;
498  fF0IsComputed = false;
499 
500  return 0;
501 }
502 
503 
504 template<class TVar, class TSideMatrix>
505 int TPZMatRed<TVar, TSideMatrix>::Zero(){
506  if(fK00) fK00->Zero();
507  fIsReduced = false;
508  fK01IsComputed = false;
509  fF0IsComputed = false;
510  fK01.Zero();
511  fK10.Zero();
512  fK11.Zero();
513  fF0.Zero();
514  fF1.Zero();
515  TPZMatrix<TVar>::Redim(fDim0+fDim1,fDim0+fDim1);
516  return 0;
517 }
518 
519 
520 template<class TVar, class TSideMatrix>
521 void TPZMatRed<TVar, TSideMatrix>::MultAdd(const TPZFMatrix<TVar> &x,
522  const TPZFMatrix<TVar> &y, TPZFMatrix<TVar> &z,
523  const TVar alpha,const TVar beta,
524  const int opt) const
525 {
526  // #warning Not functional yet. Still need to Identify all the variables
527  if(!fIsReduced)
528  {
529  LOGPZ_WARN(logger,"TPZMatRed not reduced, expect trouble")
530  TPZMatrix<TVar>::MultAdd(x,y,z,alpha,beta,opt);
531  return;
532  }
533 
534  this->PrepareZ(y,z,beta,opt);
535 
536  if(!opt)
537  {
538  if(fK01IsComputed)
539  {
540  DebugStop();
541  }
542 
543  TPZFMatrix<TVar> l_Res(fK01.Rows(), x.Cols(), 0);
544  fK01.Multiply(x,l_Res,0);
545  //fSolver->Solve(l_Res,l_Res);
546 #ifdef LOG4CXX
547  if(logger->isDebugEnabled())
548  {
549  std::stringstream sout;
550  l_Res.Print("Internal solution",sout);
551  LOGPZ_DEBUG(logger,sout.str())
552  }
553 #endif
554  TPZFMatrix<TVar> l_ResFinal(fK11.Rows(), x.Cols(), 0);
555  fK11.Multiply(x,l_ResFinal,0);
556 #ifdef LOG4CXX
557  if(logger->isDebugEnabled())
558  {
559  std::stringstream sout;
560  l_ResFinal.Print("Intermediate product l_ResFinal",sout);
561  LOGPZ_DEBUG(logger,sout.str())
562  }
563 #endif
564  fK10.MultAdd(l_Res,l_ResFinal,z,-alpha,alpha,opt);
565 #ifdef LOG4CXX
566  if(logger->isDebugEnabled())
567  {
568  std::stringstream sout;
569  z.Print("Final result z ",sout);
570  LOGPZ_DEBUG(logger,sout.str())
571  }
572 #endif
573  }
574  else
575  {
576  DebugStop();
577  }
578 }
579 
581 template<class TVar, class TSideMatrix>
582 void TPZMatRed<TVar, TSideMatrix>::DecomposeK00()
583 {
584  if(fK00->IsDecomposed())
585  {
586  return;
587  }
588 #ifdef PZDEBUG
589  {
590  REAL norm = Norm(*fK00);
591 
592 
593  if(norm < ZeroTolerance())
594  {
595  DebugStop();
596  }
597  }
598 #endif
599  TPZStepSolver<TVar> *stepsolve = dynamic_cast<TPZStepSolver<TVar> *>(fSolver.operator->());
600  TPZStepSolver<TVar> *directsolve(0);
601  if(!stepsolve)
602  {
603  DebugStop();
604  }
605  if(stepsolve->Solver() == TPZMatrixSolver<TVar>::EDirect)
606  {
607  directsolve = stepsolve;
608  }
609  if(!directsolve)
610  {
611  TPZSolver<TVar> *presolve = stepsolve->PreConditioner();
612  TPZStepSolver<TVar> *prestep = dynamic_cast<TPZStepSolver<TVar> *>(presolve);
613  if(prestep->Solver() == TPZMatrixSolver<TVar>::EDirect)
614  {
615  prestep->UpdateFrom(stepsolve->Matrix());
616  directsolve = prestep;
617  }
618  }
619  if (directsolve)
620  {
621  directsolve->Decompose();
622  std::list<int64_t> &singular = directsolve->Singular();
623  std::list<int64_t>::iterator it;
624  int nsing = singular.size();
625  if(nsing > fMaxRigidBodyModes-fNumberRigidBodyModes)
626  {
627  std::cout << "Number of rigid body modes larger than provision\n";
628  std::cout << "Number of singular modes " << nsing << std::endl;
629  std::cout << "Number of rigid body modes reserved " << fMaxRigidBodyModes << std::endl;
630  std::cout << "Rigid body modes ";
631  for (it=singular.begin(); it != singular.end(); it++) {
632  std::cout << " " << *it;
633  }
634  std::cout << std::endl;
635  //DebugStop();
636  }
637  for (it=singular.begin(); it != singular.end(); it++) {
638  if(fNumberRigidBodyModes < fMaxRigidBodyModes)
639  {
640  fK01(*it,fDim1-fMaxRigidBodyModes+fNumberRigidBodyModes) = -1.;
641  fK10(fDim1-fMaxRigidBodyModes+fNumberRigidBodyModes,*it) = -1.;
642  fK11(fDim1-fMaxRigidBodyModes+fNumberRigidBodyModes,fDim1-fMaxRigidBodyModes+fNumberRigidBodyModes) = 1.;
643  if(stepsolve != directsolve)
644  {
645  TVar diag = stepsolve->Matrix()->GetVal(*it, *it)+ (TVar)1.;
646  stepsolve->Matrix()->PutVal(*it, *it, diag);
647  }
648  }
649  fNumberRigidBodyModes++;
650  }
651  }
652  else
653  {
654  DebugStop();
655  }
656 }
657 
658 template<class TVar, class TSideMatrix>
659 void TPZMatRed<TVar, TSideMatrix>::Write(TPZStream &buf, int withclassid) const {
660  TPZMatrix<TVar>::Write(buf, withclassid);
661  {//Ints
662  buf.Write(&this->fDim0, 1);
663  buf.Write(&this->fDim1, 1);
664  }
665  {//chars
666  buf.Write(this->fIsReduced);
667  buf.Write(this->fK01IsComputed);
668  buf.Write(&this->fMaxRigidBodyModes, 1);
669  buf.Write(&this->fNumberRigidBodyModes, 1);
670  }
671  {//Aggregates
672  this->fF0.Write(buf, 0);
673  this->fF1.Write(buf, 0);
674  TPZPersistenceManager::WritePointer(this->fK00.operator ->(), &buf);
675  this->fK01.Write(buf, 0);
676  this->fK10.Write(buf, 0);
677  this->fK11.Write(buf, 0);
678  if (fSolver) {
679  if (fSolver->Matrix() != fK00) {
680  std::cout << "Error\n";
681  } else {
682  TPZPersistenceManager::WritePointer(fSolver.operator ->(), &buf);
683  //TODO Enviar o solver, atenção com a Matrix do Solver;
684  }
685 
686  } else {
687  int flag = -1;
688  buf.Write(&flag, 1);
689  }
690 
691  }
692 
693 }
694 
695 template<class TVar, class TSideMatrix>
696 void TPZMatRed<TVar, TSideMatrix>::Read(TPZStream &buf, void *context) {
697  TPZMatrix<TVar>::Read(buf, context);
698  {//Ints
699  buf.Read(&this->fDim0, 1);
700  buf.Read(&this->fDim1, 1);
701  }
702  {//chars
703  buf.Read(this->fIsReduced);
704  buf.Read(this->fK01IsComputed);
705  buf.Read(&this->fMaxRigidBodyModes, 1);
706  buf.Read(&this->fNumberRigidBodyModes, 1);
707  }
708  {//Aggregates
709  this->fF0.Read(buf, 0);
710  this->fF1.Read(buf, 0);
711  TPZAutoPointer<TPZSavable> sav = TPZPersistenceManager::GetAutoPointer(&buf);
712  TPZAutoPointer<TPZMatrix<TVar>> mat = TPZAutoPointerDynamicCast<TPZMatrix<TVar>>(sav);
713  if (sav && !mat) {
714  DebugStop();
715  }
716  fK00 = mat;
717  this->fK01.Read(buf, 0);
718  this->fK10.Read(buf, 0);
719  this->fK11.Read(buf, 0);
720  sav = TPZPersistenceManager::GetAutoPointer(&buf);
721  TPZAutoPointer<TPZMatrixSolver<TVar>>matsolv = TPZAutoPointerDynamicCast<TPZMatrixSolver<TVar>> (sav);
722  if (sav && !matsolv) {
723  DebugStop();
724  }
725  if (matsolv) {
726  fSolver = matsolv;
727  }
728  }
729 }
730 
731 #include "tpzverysparsematrix.h"
732 
733 template class TPZMatRed<float, TPZVerySparseMatrix<float> >;
734 template class TPZMatRed<float, TPZFMatrix<float> >;
735 
736 template class TPZMatRed<double, TPZVerySparseMatrix<double> >;
737 template class TPZMatRed<double, TPZFMatrix<double> >;
738 
739 template class TPZMatRed<long double, TPZVerySparseMatrix<long double> >;
740 template class TPZMatRed<long double, TPZFMatrix<long double> >;
741 
742 template class TPZMatRed<std::complex<double>, TPZVerySparseMatrix<std::complex<double> > >;
743 
744 template class TPZMatRed<std::complex<float>, TPZFMatrix<std::complex<float> > >;
745 template class TPZMatRed<std::complex<double>, TPZFMatrix<std::complex<double> > >;
746 template class TPZMatRed<std::complex<long double>, TPZFMatrix<std::complex<long double> > >;
747 
748 
749 #ifndef BORLAND
750 template class TPZRestoreClass<TPZMatRed<REAL,TPZVerySparseMatrix<REAL> > >;
751 template class TPZRestoreClass<TPZMatRed<REAL, TPZFMatrix<REAL> > >;
752 #endif
TPZMatRed::ClassId
int ClassId() const override
Saveable methods.
Definition: pzmatred.h:272
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZMatrix::Write
void Write(TPZStream &buf, int withclassid) const override
Packs the object structure in a stream of bytes.
Definition: pzmatrix.cpp:1420
TPZMatrix::Error
static int Error(const char *msg, const char *msg2=0)
Returns error messages.
Definition: pzmatrix.cpp:1402
TPZMatRed::SimetrizeMatRed
void SimetrizeMatRed()
If fK00 is simetric, only part of the matrix is accessible to external objects.
Definition: pzmatred.cpp:67
TPZMatRed::F1Red
void F1Red(TPZFMatrix< TVar > &F1)
Computes the reduced version of the right hand side .
Definition: pzmatred.cpp:171
TPZMatRed::Swap
static void Swap(int64_t *row, int64_t *col)
Swaps the row and column index.
Definition: pzmatred.h:280
TPZMatRed::UGlobal
void UGlobal(const TPZFMatrix< TVar > &U1, TPZFMatrix< TVar > &result)
Computes the complete vector based on the solution u1.
Definition: pzmatred.cpp:342
EInputFormat
Definition: pzmatrix.h:55
TPZMatRed::K11
TPZFMatrix< TVar > & K11()
Definition: pzmatred.h:125
TPZSolver
Defines a abstract class of solvers which will be used by matrix classes. Solver. ...
Definition: pzmatrix.h:32
TPZStepSolver
Defines step solvers class. Solver.
Definition: pzmganalysis.h:17
MatrixOutputFormat
MatrixOutputFormat
Defines output format.
Definition: pzmatrix.h:55
TPZMatrix::fRow
int64_t fRow
Number of rows in matrix.
Definition: pzmatrix.h:779
TPZMatRed::SetF
void SetF(const TPZFMatrix< TVar > &F)
Copies the F vector in the internal data structure.
Definition: pzmatred.cpp:140
TPZMatRed::UGlobal2
void UGlobal2(TPZFMatrix< TVar > &U1, TPZFMatrix< TVar > &result)
Definition: pzmatred.cpp:397
TPZMatrixSolver
Defines a class of matrix solvers. Solver.
Definition: pzanalysis.h:24
TPZMatrix::PrepareZ
void PrepareZ(const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar beta, const int opt) const
Is an auxiliar method used by MultiplyAdd.
Definition: pzmatrix.cpp:135
LOGPZ_WARN
#define LOGPZ_WARN(A, B)
Define log for warnings.
Definition: pzlog.h:91
TPZStepSolver::Solver
virtual TPZMatrixSolver< TVar >::MSolver Solver() override
Definition: pzstepsolver.h:65
pzmatred.h
Contains TPZMatRed class which implements a simple substructuring of a linear system of equations...
TPZMatrix::MultAdd
virtual void MultAdd(const TPZFMatrix< TVar > &x, const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar alpha=1., const TVar beta=0., const int opt=0) const
It computes z = beta * y + alpha * opt(this)*x but z and x can not overlap in memory.
Definition: pzmatrix.cpp:166
TPZMatRed::SetK00
void SetK00(TPZAutoPointer< TPZMatrix< TVar > > K00)
Sets the matrix pointer of the upper left matrix to K00.
Definition: pzmatred.cpp:134
TPZMatRed::MultAdd
void MultAdd(const TPZFMatrix< TVar > &x, const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar alpha, const TVar beta, const int opt) const override
It computes z = beta * y + alpha * opt(this)*x but z and x can not overlap in memory.
Definition: pzmatred.cpp:521
TPZMatRed::fK01IsComputed
bool fK01IsComputed
Is true if has been computed and overwritten .
Definition: pzmatred.h:259
TPZMatRed::SetSolver
void SetSolver(TPZAutoPointer< TPZMatrixSolver< TVar > > solver)
Definition: pzmatred.cpp:125
TPZMatRed::s
virtual TVar & s(const int64_t row, const int64_t col) override
The operators check on the bounds if the DEBUG variable is defined.
Definition: pzmatred.cpp:113
TPZMatRed::Print
void Print(const char *name=NULL, std::ostream &out=std::cout, const MatrixOutputFormat=EFormatted) const override
Prints the object data structure.
Definition: pzmatred.cpp:454
TPZPersistenceManager::GetAutoPointer
static TPZAutoPointer< TPZSavable > GetAutoPointer(const int64_t &objId)
Definition: TPZPersistenceManager.cpp:417
TPZMatRed::fK10
TSideMatrix fK10
Definition: pzmatred.h:247
TPZStepSolver::Singular
std::list< int64_t > & Singular()
returns the equations for which the equations had zero pivot
Definition: pzstepsolver.h:71
TPZMatRed::Read
void Read(TPZStream &buf, void *context) override
read objects from the stream
Definition: pzmatred.cpp:696
TPZMatRed::fF0
TPZFMatrix< TVar > fF0
Right hand side or force matrix.
Definition: pzmatred.h:250
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZStepSolver::Decompose
virtual void Decompose() override
Decompose the system of equations if a direct solver is used.
Definition: pzstepsolver.cpp:57
TPZStepSolver::UpdateFrom
virtual void UpdateFrom(TPZAutoPointer< TPZMatrix< TVar > > matrix) override
Updates the values of the current matrix based on the values of the matrix.
Definition: pzstepsolver.h:81
Norm
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
TPZMatRed::fSolver
TPZAutoPointer< TPZMatrixSolver< TVar > > fSolver
Solution method for inverting .
Definition: pzmatred.h:242
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZVerySparseMatrix
Implements a matrix whose nonzero elements are stored in binary tree. Matrix.
Definition: pzfmatrix.h:33
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
ELU
Definition: pzmatrix.h:52
TPZMatRed::fF0IsComputed
bool fF0IsComputed
Is true if has been computed and overwritten .
Definition: pzmatred.h:262
TPZMatRed::fIsReduced
bool fIsReduced
Is true if the declared dimension of the matrix is fDim0.
Definition: pzmatred.h:256
TPZMatRed::fK00
TPZAutoPointer< TPZMatrix< TVar > > fK00
Stiffnes matrix.
Definition: pzmatred.h:239
TPZMatRed::fF1
TPZFMatrix< TVar > fF1
Definition: pzmatred.h:250
TPZMatrix::Read
void Read(TPZStream &buf, void *context) override
Unpacks the object structure from a stream of bytes.
Definition: pzmatrix.cpp:1413
TPZMatRed::fK11
TPZFMatrix< TVar > fK11
Full Stiffnes matrix.
Definition: pzmatred.h:245
TPZMatrix::Redim
virtual int Redim(const int64_t newRows, const int64_t newCols)
Redimensions the matrix reinitializing it with zero.
Definition: pzmatrix.h:289
TPZFMatrix
Full matrix class. Matrix.
Definition: pzfmatrix.h:32
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZMatRed::Redim
int Redim(const int64_t dim, const int64_t dim00) override
Redim: Set the dimension of the complete matrix and reduced matrix.
Definition: pzmatred.cpp:481
TPZMatrix::Dim
virtual int64_t Dim() const
Returns the dimension of the matrix if the matrix is square.
Definition: pzmatrix.h:892
TPZMatRed::DecomposeK00
void DecomposeK00()
Decompose K00 and adjust K01 and K10 to reflect rigid body modes.
Definition: pzmatred.cpp:582
TPZMatrix::fCol
int64_t fCol
Number of cols in matrix.
Definition: pzmatrix.h:781
TPZMatRed::PutVal
virtual int PutVal(const int64_t row, const int64_t col, const TVar &value) override
Put and Get values without bounds checking these methods are faster than "Put" e "Get" if DEBUG is de...
Definition: pzmatred.cpp:88
TPZMatRed::fNumberRigidBodyModes
int fNumberRigidBodyModes
Number of rigid body modes identified during the decomposition of fK00.
Definition: pzmatred.h:268
TPZMatRed
Implements a simple substructuring of a linear system of equations, composed of 4 submatrices...
Definition: pzmatred.h:34
TPZStepSolver::PreConditioner
TPZSolver< TVar > * PreConditioner()
access method to the preconditioner
Definition: pzstepsolver.h:104
TPZMatRed::K00
TPZAutoPointer< TPZMatrix< TVar > > K00()
Definition: pzmatred.h:112
TPZMatRed::Write
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: pzmatred.cpp:659
TPZMatRed::fK01
TSideMatrix fK01
Definition: pzmatred.h:247
TPZMatRed::K11Reduced
void K11Reduced(TPZFMatrix< TVar > &K11, TPZFMatrix< TVar > &F1)
Computes the K11 reduced .
Definition: pzmatred.cpp:223
TPZMatrixSolver::Matrix
TPZAutoPointer< TPZMatrix< TVar > > Matrix() const
Returns a pointer to TPZMatrix<>
Definition: pzsolve.h:138
TPZPersistenceManager.h
TPZMatRed::IsSimetric
virtual int IsSimetric() const override
returns 1 or 0 depending on whether the fK00 matrix is zero or not
Definition: pzmatred.cpp:61
EMathematicaInput
Definition: pzmatrix.h:55
pzstepsolver.h
Contains TPZStepSolver class which defines step solvers class.
TPZMatRed::fDim1
int64_t fDim1
Definition: pzmatred.h:253
tpzverysparsematrix.h
Contains TPZVerySparseMatrix class which implements a matrix whose nonzero elements are stored in bin...
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZMatRed::fDim0
int64_t fDim0
Stores matricess dimensions.
Definition: pzmatred.h:253
TPZMatRed::F1
TPZFMatrix< TVar > & F1()
Definition: pzmatred.h:130
TPZMatRed::GetVal
virtual const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzmatred.cpp:101
TPZMatRed::fMaxRigidBodyModes
int fMaxRigidBodyModes
Number of rigid body modes foreseen in the computational mesh.
Definition: pzmatred.h:265
TPZMatRed::~TPZMatRed
~TPZMatRed()
Simple destructor.
Definition: pzmatred.cpp:57
ZeroTolerance
REAL ZeroTolerance()
Returns the tolerance to Zero value. Actually: .
Definition: pzreal.h:633
TPZFMatrix::PutVal
int PutVal(const int64_t row, const int64_t col, const TVar &value) override
Put values without bounds checking This method is faster than "Put" if DEBUG is defined.
Definition: pzfmatrix.h:548
TPZFMatrix::GetVal
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
TPZPersistenceManager::WritePointer
static void WritePointer(const TPZSavable *obj, TPZStream *stream)
Definition: TPZPersistenceManager.cpp:182
TPZMatRed::Zero
virtual int Zero() override
This method will zero all submatrices associated with this reducable matrix class.
Definition: pzmatred.cpp:505
TPZRestoreClass
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
TPZMatrix::SolveDirect
virtual int SolveDirect(TPZFMatrix< TVar > &F, const DecomposeType dt, std::list< int64_t > &singular)
Solves the linear system using Direct methods.
Definition: pzmatrix.cpp:710
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
TPZMatRed::TPZMatRed
TPZMatRed()
Simple constructor.
Definition: pzmatred.cpp:31
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
TPZMatRed::U1
void U1(TPZFMatrix< TVar > &F)
Returns the second vector, inverting K11.
Definition: pzmatred.cpp:289
TPZMatrix
Root matrix class (abstract). Matrix.
Definition: pzmatrix.h:60
TPZAutoPointer
This class implements a reference counter mechanism to administer a dynamically allocated object...
Definition: TPZPersistenceManager.h:14
Generated on Sun Aug 16 2020 11:17:05 for NeoPZ by   doxygen 1.8.13