Commits

Lisandro Dalcin committed f886724

Fix Neumann demo to properly handle the constant null space

Comments (0)

Files changed (1)

 
 PetscReal Exact(PetscReal x[3])
 {
-  return sin(pi*x[0]) + sin(pi*x[1]) + sin(pi*x[2]);
+  return sin(2*pi*x[0]) + sin(2*pi*x[1]) + sin(2*pi*x[2]);
 }
 
 PetscReal Forcing(PetscReal x[3])
 {
-  return pi*pi * Exact(x);
+  return 4*pi*pi * Exact(x);
 }
 
 PetscReal Flux(PetscInt dir,PetscInt side)
 {
-  return -pi;
+  return (side?+1:-1) * 2*pi;
 }
 
 PETSC_STATIC_INLINE
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "MassVector"
+PetscErrorCode MassVector(IGAPoint p,PetscScalar *V,void *ctx)
+{
+  PetscInt a,nen = p->nen;
+  const PetscReal *N0 = (typeof(N0)) p->shape[0];
+  for (a=0; a<nen; a++) V[a] = N0[a];
+  return 0;
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "Error"
 PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
 {
   IGAPointFormPoint(p,x);
   PetscScalar u;
   IGAPointFormValue(p,U,&u);
-  PetscReal e = PetscAbsScalar(u) -  Exact(x);
+  PetscScalar e = u -  Exact(x);
   S[0] = e*e;
   return 0;
 }
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
-  {
-    PetscReal vmin; /* this is a hack */
-    ierr = VecMin(x,0,&vmin);CHKERRQ(ierr);
-    ierr = VecShift(x,-vmin);CHKERRQ(ierr);
-  }
+
+  Vec Q;
+  ierr = IGACreateVec(iga,&Q);CHKERRQ(ierr);
+  ierr = IGASetFormVector(iga,MassVector,NULL);CHKERRQ(ierr);
+  ierr = IGAComputeVector(iga,Q);CHKERRQ(ierr);
+  PetscScalar mean;
+  ierr = VecDot(Q,x,&mean);CHKERRQ(ierr);
+  ierr = VecShift(x,-mean);CHKERRQ(ierr);
 
   PetscScalar error;
   ierr = IGAComputeScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.