Επιστροφή στην αρχική σελίδα...

O Αλγόριθμος Επίλυσης Skyline

Κατά την επίλυση με την Μέθοδο Πεπερασμένων Στοιχείων, είναι απαραίτητη η επίλυση μέσω της τριγωνοποίησης του συστήματος στιβαρότητας:

[K].u = F

Όπου [K] το μητρώο στιβαρότητας, u το (άγνωστο) διάνυσμα-στήλη μετατοπίσεων και F το διάνυσμα-στήλη των ασκούμενων φορτίων.

Θα ήταν επιθυμητός ο σχηματισμός ενός μητρώου [N] όπου η πίσω-αντικατάσταση θα γινόταν με την χρήση αυτού του μητρώου. Οι αλγόριθμοι που βασίζονται στην απαλοιφή Gauss είναι ακατάλληλοι καθώς κατά την διαδικασία της τριγωνοποίησης «πειράζεται» και το διάνυσμα των σταθερών όρων (F στην περίπτωσή μας). Λύση δίνουν οι μέθοδοι LU decomposition (βλέπε Chapra & Canale, Numerical Methods for Engineers, McGraw Hill). Το πλεονέκτημα της LU decomp. είναι πως γίνεται μία φορά η επίλυση και σε κάθε iteration ή loading step γίνεται απλά μία πίσω αντικατάσταση. Η οικονομία σε υπολογισμούς είναι σημαντική έως θαυματουργή!

Επιπλέον τα μητρώα στιβαρότητας είναι συμμετρικά στις περισσότερες περιπτώσεις (μία εξαίρεση είναι η χρήση μη-γραμμικού μητρώου στιβαρότητας με μη-συσχετισμένο νόμο ροής). Ακόμα, με διαδικασία βελτιστοποίησης κατά την διαμέριση σε πεπερασμένα στοιχεία/ κόμβους, το μητρώο είναι συγκεντρωμένο γύρω από την κύρια διαγώνιο με ένα ημιεύρος HBW. (αν έχουμε π.χ. 3 βαθμούς ελευθερίας ανά κόμβο και η μέγιστη διαφορά μεταξύ αρίθμησης κόμβων που ενώνει ένα στοιχείο είναι n, το HBW θα είναι HBW=3n)

Τέλος είναι δυνατή η εκμετάλλευση του «προφίλ» του μητρώου που θυμίζει γραμμή από ουρανοξύστες (skyline), λόγω των πολλών μηδενικών στοιχείων.

Αν είναι γνωστό το HBW και έχει σχηματιστεί ένα διάνυσμα που περιέχει την Skyline ως μέγιστο μη-μηδενικό στοιχείο κατά στήλες, η επίλυση γίνεται με τον παρακάτω κώδικα (απόσπασμα από το CEFEAS):

/***********************************/
/** Skyline Solving               **/
/***********************************/
void skysolve(int n,int mode)
{
 int i,j,k,r;
 progress progress;

/***************************************/
/** Elimination by LDL` Factorisation **/
/***************************************/
 if(mode==0){
  progress.open("Forward Elimination",n);
  for(j=2;j<=n;j++){
   progress.refresh(j);
   for(i=Skyline[j]+1;i<j;i++){
    for(r=max(Skyline[i],Skyline[j]);r<i;r++){
     K[pos(i,j-i+1)]-=K[pos(r,i-r+1)]*K[pos(r,j-r+1)]; //Formulate gij
    }
   }
   for(i=Skyline[j];i<j;i++){
    U[i]=K[pos(i,j-i+1)];
    if(K[pos(i,1)]!=0.000)  
     K[pos(i,j-i+1)]/=K[pos(i,1)]; //Formulate lij
    else{
     PRINT("Forward Elimination: Division by Zero \n");
    }
   }
   for(r=Skyline[j];r<j;r++)
    K[pos(j,1)]-=K[pos(r,j-r+1)]*U[r]; //Formulate djj 
  } 
  progress.close(); 
 }

/*************************************/
/** Back Substitution Process       **/
/*************************************/
 progress.open("Back Substitution",n);
 for(i=2;i<=n;i++){
  progress.refresh(i);
  for(r=Skyline[i];r<i;r++)
   Q[i]-=K[pos(r,i-r+1)]*Q[r]; 
 } 
 for(i=1;i<=n;i++){
  if(K[pos(i,1)]!=0.000) 
   Q[i]/=K[pos(i,1)];
   else{
    PRINT("Back Substitution: Division by Zero \n");
   }
 }
 for(i=n;i>=2;i--){
  progress.refresh(i);
  for(r=Skyline[i];r<i;r++){
   Q[r]=Q[r]-K[pos(r,i-r+1)]*Q[i];
  }            
 } 
 progress.close();
 for(i=1;i<=n;i++)
  U[i]=Q[i];
}                 

 Αγνοείστε τα PRINT, progress. Το pos(i,j) απεικονίζει την γραμμή i, στήλη j, σε ένα μονοδιάστατο αριθμό για πρόσβαση στο μονοδιάστατο μητρώο K. Q=F. Το Skyline[] μετράει από την γραμμή 1 και όχι από την διαγώνιο. Δηλαδή αν Skyline[5]=1 σημαίνει πως η 5η γραμμή φτάνει ως πάνω-πάνω, ενώ Skyline[6]=6 σημαίνει πως στην έκτη στήλη έχουμε μόνο στοιχείο στην διαγώνιο. Οι πίνακες δεν είναι zero-based, ξεκινάμε από το 1.

Επιστροφή στην αρχική σελίδα...