Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Neuman Boundary Conditions

2011-07-21
2013-03-14
  • Thomas Plehn
    Thomas Plehn
    2011-07-21

    I am new to ofeli, and I am very curious about the appearantly not working neuman condition, which I prescribed for specific sides (lines) using the code in the domain.addLine(…,dc,nc=5). My Prescription file has the following line in it:
    <BoundaryForce code="5" dof="1">100<BoundaryForce>
    or do I do something wrong?

    to exclude the fact, that it has something to do with the solver, the solver code is attached:

    if (output_flag > 1)
          cout << "Reading mesh data …\n";
       Mesh ms("mymesh.m");

       saveGmsh("test.gmsh",ms);

       Prescription p(ms,"truss.pre");
       if (output_flag > 1)
          cout << ms;

    // Declare problem data (matrix, rhs, boundary conditions, body forces)
       if (output_flag > 1)
          cout << "Allocating memory for matrix and R.H.S. …\n";
       SkSMatrix<double> A(ms);
       Vect<double> b(ms.getNbDOF());

    // Read boundary conditions, body and boundary forces
       if (output_flag > 1)
         cout << "Reading boundary conditions …" << endl;
       NodeVect<double> bc(ms,ms.getNbDOF());
       p.get(BOUNDARY_CONDITION,bc);
       if (output_flag > 1)
         cout << "Reading body forces …" << endl;
       NodeVect<double> body_f(ms,ms.getNbDOF());
       p.get(BODY_FORCE,body_f);
       if (output_flag > 1)
         cout << "Reading Boundary Tractions …" << endl;
       SideVect<double> bound_f(ms,2);
       p.get(BOUNDARY_FORCE,bound_f);

    // Loop over elements
    // ------------

       if (output_flag > 1)
          cout << "Looping over elements …\n";

       MeshElements(ms) {
          Elas2DT3 eq(theElement);
      eq.PlaneStress();
          eq.Deviator();
          eq.Dilatation();
          eq.BodyRHS(body_f);
          eq.ElementAssembly(A);
          eq.ElementAssembly(b);
       }

    // Loop over sides
    // ----------

       if (output_flag > 1)
         cout << "Looping over sides …\n";
       MeshSides(ms) {
          Elas2DT3 eq(theSide);
      eq.PlaneStress();
          eq.BoundaryRHS(bound_f);
          eq.SideAssembly(b);
      eq.SideAssembly(A);
       }

    // Take account for boundary conditions and solve system
    // -----------------------------------

       if (output_flag > 1)
         cout << "Imposing boundary conditions …\n";
       A.Prescribe(ms,b,bc);
       A.Factor();
       A.Solve(b);