diff --git a/CHANGELOG.txt b/CHANGELOG.txt index d0d48141..d7ba324b 100644 --- a/CHANGELOG.txt +++ b/CHANGELOG.txt @@ -17,7 +17,7 @@ TODO for 0.5.0 - "Split package, multiple jars provide the same package" - [WARNING] bootstrap class path not set in conjunction with -source 8 when using: mvn clean install -P on-jdk-9-plus -- CI: Fail build on warnbing?! +- CI: Fail build on warning?! - CI: remove "temurin"...?! @@ -42,6 +42,10 @@ ODE16.2: Much better. Floating card is actually fine (in theory). But3rd row col 0.5.0 (unreleased) ===== +- Port updates until 0.16.2. + This includes some libCCD updates missing fro the 0.15.1 update. + This excludes improved solution finders with anything from *Cooperative* and ThreadedEquationSolverLDLT*. + [#97](https://github.com/tzaeschke/ode4j/pull/97) - Added GutHub Actions CI builds for Java 8 and 9. [#95](https://github.com/tzaeschke/ode4j/pull/95) - Added default logger for demos and tests. [#94](https://github.com/tzaeschke/ode4j/pull/94) - Added Android API compliance checker, now for API level 24 (the lowest that ode4j passed without changes) @@ -55,7 +59,7 @@ ODE16.2: Much better. Floating card is actually fine (in theory). But3rd row col (Cleanup, part of this was done with the move to 0.15.1) - Update to ODE 0.15.1 [#86](https://github.com/tzaeschke/ode4j/pull/86) - This includes increased stability, e.g. the DemoCards won't collapse anymore when 30 levels high. - - **** WARNING **** : Still missing: libccd and trimesh collider updates. + - **** WARNING **** : Still missing: trimesh collider updates. - Changed CI build requirements: - Removed openjdk7: Not required. source compatibility is still set to 7, but should be safe to change to 8 now (Android Nougat is mostly 8 compatible), this covers > 80% users. diff --git a/README.md b/README.md index 6299cd47..2f7307c0 100644 --- a/README.md +++ b/README.md @@ -68,16 +68,20 @@ News * This is the last release built with Java 6. +Basic Usage Tips +================ + * Use `World.quickStep(...)` instead of `World.step()`. The latter is slower and appears to be less stable. + * Make sure to set `Common.dNODEBUG = true` for best performance. + * Avoid using `core-cpp` and ignore demos in `demo-cpp`. + + Legal ===== ode4j: -Copyright (c) 2009-2017 Tilmann Zäschke . +Copyright (c) 2009-2023 Tilmann Zäschke . All rights reserved. - - - Like the original ODE, ode4j is licensed under LGPL v2.1 and BSD 3-clause. Choose whichever license suits your needs. diff --git a/TODO.txt b/TODO.txt index 44545f6d..ac103e67 100644 --- a/TODO.txt +++ b/TODO.txt @@ -5,13 +5,24 @@ mvn versions:display-dependency-updates mvn versions:display-plugin-updates +Things left out for 0.5.0 (missing from 0.15.1 vs 0.16.2) +========================================================= + Not ported: + - resource_control.* + - simple_cooperation.* + - Updates to DxWorldProcessContext + - ThreadedEquationSolverLDLT* and *Cooperative* is missing from: + - fastldltfactor(_impl).* + - fastlsolve(_impl).* + - fastltsolve(_impl).*: + - fastvecscale(_impl).*: + - objects.h/.cpp / obstack.h/.cpp + Java Version ============ - We should aim for Oreo or later (Android 8.0, released 2017, API Level 26: https://apilevels.com/) -- Java 8 should be fine. -- Try animal-sniffer plugin for higher versions, e.g. Java 9 - More accurate (respects desugaring by D8): https://github.com/open-toast/gummy-bears +- We use animal-sniffer (respects desugaring by D8): https://github.com/open-toast/gummy-bears - Android Runtime replaced Dalvik as default Runtime in Android 5: https://en.wikipedia.org/wiki/Android_version_history - Android 8 (Oreo) is about 8% share: @@ -24,9 +35,6 @@ Java Version - Android 11 support "up to" Java 13 (Question: complete 13 or partial 13??) - 2018: https://jakewharton.com/androids-java-9-10-11-and-12-support/ -- Java 8 support: https://developer.android.com/studio/write/java8-support-table -- Java 8 support details: https://developer.android.com/studio/write/java8-support - MArch 2023: General cleanup: @@ -35,8 +43,12 @@ General cleanup: - Move tests to core/test - update DVector3 safeNormalize to use max() and then compare to 0. ALso see: https://github.com/JOML-CI/JOML/issues/66 - replace DVector3.clone() with dVector3.re() or copy() +- Create separate distribution (maven jar) with modules? -> Check how modul-projects have trouble with including + non-module projects... + -> https://www.baeldung.com/java-9-modularity + -> https://stackoverflow.com/questions/40490520/what-do-i-need-to-build-jdk-9-project-with-non-modular-dependencies-using-maven - reformat everything? At least indentation? -- CI builds +- Consider porting stuff in "contrib" July 6th 2009: @@ -49,7 +61,7 @@ Some minor things I found in the code: - http://opende.sourceforge.net/wiki/index.php/Joint_Param_Documentation This indicate that PU & PR joint have not parameters, I don't think that's correct? - GeomTransform is obsolete since 0.6 (according to the Wiki), but it is still used in several demos: - deom_boxstack, demo_heightfield, demo_jointPU, demo_motion, demo_spacestress. + demo_boxstack, demo_heightfield, demo_jointPU, demo_motion, demo_spacestress. I haven't used it myself, so is it still required? diff --git a/core-cpp/src/main/java/org/ode4j/cpp/internal/ApiCppMathMatrix.java b/core-cpp/src/main/java/org/ode4j/cpp/internal/ApiCppMathMatrix.java index 819e875f..1e8d94d4 100644 --- a/core-cpp/src/main/java/org/ode4j/cpp/internal/ApiCppMathMatrix.java +++ b/core-cpp/src/main/java/org/ode4j/cpp/internal/ApiCppMathMatrix.java @@ -218,11 +218,21 @@ void dSolveL1T (final DMatrix3C L, DVector3 b, int n, int nskip) { throw new UnsupportedOperationException(); } + /* in matlab syntax: a(1:n) = a(1:n) .* d(1:n) + */ - /** in matlab syntax: a(1:n) = a(1:n) .* d(1:n) */ + // ODE_API void dScaleVector (dReal *a, const dReal *d, int n); + void dScaleVector (DVector3 a, final DVector3C d, int n) { + throw new UnsupportedOperationException(); + } - //ODE_API - // void dVectorScale (double *a, final double *d, int n) { + /* The function is an alias for @c dScaleVector. + * It has been deprecated because of a wrong naming schema used. + */ + //ODE_API_DEPRECATED + // ODE_API + // void dVectorScale (dReal *a, const dReal *d, int n) + @Deprecated void dVectorScale (DVector3C a, final DVector3 d, int n) { throw new UnsupportedOperationException(); } diff --git a/core/src/main/java/org/ode4j/ode/DMatrix.java b/core/src/main/java/org/ode4j/ode/DMatrix.java index 64eb90f8..bdc750f0 100644 --- a/core/src/main/java/org/ode4j/ode/DMatrix.java +++ b/core/src/main/java/org/ode4j/ode/DMatrix.java @@ -318,11 +318,26 @@ public static void dFactorLDLT (double[] A, double[] d, int n, int nskip) { // public void dSolveL1T (const dReal *L, dReal *b, int n, int nskip); + /* in matlab syntax: a(1:n) = a(1:n) .* d(1:n) + */ + + /** + * In matlab syntax: a(1:n) = a(1:n) .* d(1:n) + * @param a a + * @param d d + */ + // ODE_API + public static void dScaleVector (DVector3 a, DVector3C d) { + a.scale(d); + } + /** - * In matlab syntax: a(1:n) = a(1:n) .* d(1:n) + * The function is an alias for @c dScaleVector. + * It has been deprecated because of a wrong naming schema used. * @param a a * @param d d */ + @Deprecated // deprecated in ODE public static void dVectorScale (DVector3 a, DVector3C d) { a.scale(d); } diff --git a/core/src/main/java/org/ode4j/ode/DWorld.java b/core/src/main/java/org/ode4j/ode/DWorld.java index 24d5e0f3..510c223c 100644 --- a/core/src/main/java/org/ode4j/ode/DWorld.java +++ b/core/src/main/java/org/ode4j/ode/DWorld.java @@ -185,7 +185,9 @@ public interface DWorld { */ double getCFM() ; - public static int dWORLDSTEP_THREADCOUNT_UNLIMITED = 0; + // TZ: This should be in threading.h + static int dTHREADING_THREAD_COUNT_UNLIMITED = 0; + public static int dWORLDSTEP_THREADCOUNT_UNLIMITED = dTHREADING_THREAD_COUNT_UNLIMITED; /** * Set maximum threads to be used for island stepping diff --git a/core/src/main/java/org/ode4j/ode/internal/CollideTrimeshBoxOld.java b/core/src/main/java/org/ode4j/ode/internal/CollideTrimeshBoxOld.java deleted file mode 100644 index db25acc1..00000000 --- a/core/src/main/java/org/ode4j/ode/internal/CollideTrimeshBoxOld.java +++ /dev/null @@ -1,1646 +0,0 @@ -/************************************************************************* - * * - * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith. * - * All rights reserved. Email: russ@q12.org Web: www.q12.org * - * Open Dynamics Engine 4J, Copyright (C) 2009-2014 Tilmann Zaeschke * - * All rights reserved. Email: ode4j@gmx.de Web: www.ode4j.org * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of EITHER: * - * (1) The GNU Lesser General Public License as published by the Free * - * Software Foundation; either version 2.1 of the License, or (at * - * your option) any later version. The text of the GNU Lesser * - * General Public License is included with this library in the * - * file LICENSE.TXT. * - * (2) The BSD-style license that is included with this library in * - * the file ODE-LICENSE-BSD.TXT and ODE4J-LICENSE-BSD.TXT. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files * - * LICENSE.TXT, ODE-LICENSE-BSD.TXT and ODE4J-LICENSE-BSD.TXT for more * - * details. * - * * - *************************************************************************/ -package org.ode4j.ode.internal; - -import org.ode4j.math.DMatrix3; -import org.ode4j.math.DMatrix3C; -import org.ode4j.math.DVector3; -import org.ode4j.math.DVector3C; -import org.ode4j.math.DVector4; -import org.ode4j.math.DVector4C; -import org.ode4j.ode.DAABBC; -import org.ode4j.ode.DColliderFn; -import org.ode4j.ode.DContactGeom; -import org.ode4j.ode.DContactGeomBuffer; -import org.ode4j.ode.DGeom; -import org.ode4j.ode.internal.cpp4j.java.RefBoolean; -import org.ode4j.ode.internal.cpp4j.java.RefDouble; -import org.ode4j.ode.internal.cpp4j.java.RefInt; -import org.ode4j.ode.internal.gimpact.GimDynArrayInt; -import org.ode4j.ode.internal.gimpact.GimTrimesh; -import org.ode4j.ode.internal.gimpact.GimGeometry.aabb3f; -import org.ode4j.ode.internal.gimpact.GimGeometry.vec3f; -import org.ode4j.ode.internal.trimesh.DxTriMesh; - -import static org.ode4j.ode.OdeConstants.*; -import static org.ode4j.ode.internal.Common.*; - -/************************************************************************* - * * - * Triangle-box collider by Alen Ladavac and Vedran Klanac. * - * Ported to ODE by Oskari Nyman. * - * Ported to Java by Tilmann Zaeschke. * - * * - *************************************************************************/ -public class CollideTrimeshBoxOld implements DColliderFn { - - @Override - public int dColliderFn(DGeom o1, DGeom o2, int flags, - DContactGeomBuffer contacts) { - return dCollideBTL((DxTriMesh)o1, (DxBox)o2, flags, contacts, 1); - } - - // #if dTRIMESH_ENABLED - - - // static void - // GenerateContact(int in_Flags, dContactGeom* in_Contacts, int in_Stride, - // dxGeom* in_g1, dxGeom* in_g2, int TriIndex, - // const dVector3 in_ContactPos, const dVector3 in_Normal, dReal in_Depth, - // int& OutTriCount); - - - // largest number, double or float - // #if defined(dSINGLE) - // #define MAXVALUE FLT_MAX - // #else - // #define MAXVALUE DBL_MAX - // #endif - private static final double MAXVALUE = Double.MAX_VALUE; - - - // dVector3 - // r=a-b - // #define SUBTRACT(a,b,r) do{ \ - // (r)[0]=(a)[0] - (b)[0]; \ - // (r)[1]=(a)[1] - (b)[1]; \ - // (r)[2]=(a)[2] - (b)[2]; }while(0) - private void SUBTRACT(DVector3C a, DVector3C b, DVector3 r) { - r.eqDiff(a, b); - } - - - // dVector3 - // a=b - // #define SET(a,b) do{ \ - // (a)[0]=(b)[0]; \ - // (a)[1]=(b)[1]; \ - // (a)[2]=(b)[2]; }while(0) - private void SET(DVector3 a, DVector3C b) { - a.set(b); - } - - - // dMatrix3 - // a=b - // #define SETM(a,b) do{ \ - // (a)[0]=(b)[0]; \ - // (a)[1]=(b)[1]; \ - // (a)[2]=(b)[2]; \ - // (a)[3]=(b)[3]; \ - // (a)[4]=(b)[4]; \ - // (a)[5]=(b)[5]; \ - // (a)[6]=(b)[6]; \ - // (a)[7]=(b)[7]; \ - // (a)[8]=(b)[8]; \ - // (a)[9]=(b)[9]; \ - // (a)[10]=(b)[10]; \ - // (a)[11]=(b)[11]; }while(0) - // /** a.set(b); */ - // private void SETM(DMatrix3 a, DMatrix3C b) { - // a.set(b); - // } - - - // dVector3 - // r=a+b - // #define ADD(a,b,r) do{ \ - // (r)[0]=(a)[0] + (b)[0]; \ - // (r)[1]=(a)[1] + (b)[1]; \ - // (r)[2]=(a)[2] + (b)[2]; }while(0) - private void ADD(DVector3C a, DVector3C b, DVector3 r) { - r.eqSum(a, b); - } - - - // dMatrix3, int, dVector3 - // v=column a from m - // #define GETCOL(m,a,v) do{ \ - // (v)[0]=(m)[(a)+0]; \ - // (v)[1]=(m)[(a)+4]; \ - // (v)[2]=(m)[(a)+8]; }while(0) - private void GETCOL(DMatrix3C m, int a, DVector3 v) { - if (a==0) { v.set(m.get00(), m.get10(), m.get20()); } - else if (a==1) { v.set(m.get01(), m.get11(), m.get21()); } - else if (a==2) { v.set(m.get02(), m.get12(), m.get22()); } - else throw new IllegalArgumentException("col=" + a); - } - - - // dVector4, dVector3 - // distance between plane p and point v - // #define POINTDISTANCE(p,v) \ - // ( p[0]*v[0] + p[1]*v[1] + p[2]*v[2] + p[3] ) - private double POINTDISTANCE(DVector4C p, DVector3C v) { - return p.get0()*v.get0() + p.get1()*v.get1() + p.get2()*v.get2() + p.get3(); - } - - - // dVector4, dVector3, dReal - // construct plane from normal and d - // #define CONSTRUCTPLANE(plane,normal,d) do{ \ - // plane[0]=normal[0];\ - // plane[1]=normal[1];\ - // plane[2]=normal[2];\ - // plane[3]=d; }while(0) - private void CONSTRUCTPLANE(DVector4 plane, DVector3C n, double d) { - plane.set( n.get0(), n.get1(), n.get2(), d); - } - - - // dVector3 - // length of vector a - // #define LENGTHOF(a) \ - // dSqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]) - private double LENGTHOF(DVector3C a) { - return a.length(); - } - - - private class sTrimeshBoxColliderData - { - //sTrimeshBoxColliderData(): m_iBestAxis(0), m_iExitAxis(0), m_ctContacts(0) {} - sTrimeshBoxColliderData() { m_iBestAxis=(0); m_iExitAxis=(0); m_ctContacts=new RefInt(0); } - - // void SetupInitialContext(dxTriMesh *TriMesh, dxGeom *BoxGeom, - // int Flags, dContactGeom* Contacts, int Stride); - // int TestCollisionForSingleTriangle(int ctContacts0, int Triint, - // dVector3 dv[3], bool &bOutFinishSearching); - // - // bool _cldTestNormal(dReal fp0, dReal fR, dVector3 vNormal, int iAxis); - // bool _cldTestFace(dReal fp0, dReal fp1, dReal fp2, dReal fR, dReal fD, - // dVector3 vNormal, int iAxis); - // bool _cldTestEdge(dReal fp0, dReal fp1, dReal fR, dReal fD, - // dVector3 vNormal, int iAxis); - // bool _cldTestSeparatingAxes(const dVector3 &v0, const dVector3 &v1, const dVector3 &v2); - // void _cldClipping(const dVector3 &v0, const dVector3 &v1, const dVector3 &v2, int TriIndex); - // void _cldTestOneTriangle(const dVector3 &v0, const dVector3 &v1, const dVector3 &v2, int TriIndex); - - // box data - final DMatrix3 m_mHullBoxRot = new DMatrix3(); - final DVector3 m_vHullBoxPos = new DVector3(); - final DVector3 m_vBoxHalfSize = new DVector3(); - - // mesh data - final DVector3 m_vHullDstPos = new DVector3(); - - // global collider data - final DVector3 m_vBestNormal = new DVector3(); - //dReal m_fBestDepth; - double m_fBestDepth; - int m_iBestAxis; - @SuppressWarnings("unused") - int m_iExitAxis; - final DVector3 m_vE0 = new DVector3(), m_vE1 = new DVector3(), m_vE2 = new DVector3(), m_vN = new DVector3(); - - // global info for contact creation - int m_iFlags; - // dContactGeom *m_ContactGeoms; - DContactGeomBuffer m_ContactGeoms; - int m_iStride; - // dxGeom *m_Geom1; - // dxGeom *m_Geom2; - DxGeom m_Geom1; - DxGeom m_Geom2; - RefInt m_ctContacts; - // }; - - // Test normal of mesh face as separating axis for intersection - //boolean sTrimeshBoxColliderData::_cldTestNormal(dReal fp0, dReal fR, dVector3 vNormal, int iAxis) - boolean _cldTestNormal(double fp0, double fR, DVector3C vNormal, int iAxis) - { - // calculate overlapping interval of box and triangle - double fDepth = fR+fp0; - - // if we do not overlap - if ( fDepth<0 ) { - // do nothing - return false; - } - - // calculate normal's length - double fLength = LENGTHOF(vNormal); - // if long enough - if ( fLength > 0.0f ) { - - double fOneOverLength = 1.0f/fLength; - // normalize depth - fDepth = fDepth*fOneOverLength; - - // get minimum depth - if (fDepth < m_fBestDepth) { - // m_vBestNormal[0] = -vNormal[0]*fOneOverLength; - // m_vBestNormal[1] = -vNormal[1]*fOneOverLength; - // m_vBestNormal[2] = -vNormal[2]*fOneOverLength; - m_vBestNormal.set(vNormal).scale(-fOneOverLength); - m_iBestAxis = iAxis; - //dAASSERT(fDepth>=0); - m_fBestDepth = fDepth; - } - } - - return true; - } - - - - - // Test box axis as separating axis - // bool sTrimeshBoxColliderData::_cldTestFace(dReal fp0, dReal fp1, dReal fp2, dReal fR, dReal fD, - // dVector3 vNormal, int iAxis) - boolean _cldTestFace(double fp0, double fp1, double fp2, double fR, double fD, - DVector3 vNormal, int iAxis) - { - double fMin, fMax; - - // find min of triangle interval - if ( fp0 < fp1 ) { - if ( fp0 < fp2 ) { - fMin = fp0; - } else { - fMin = fp2; - } - } else { - if( fp1 < fp2 ) { - fMin = fp1; - } else { - fMin = fp2; - } - } - - // find max of triangle interval - if ( fp0 > fp1 ) { - if ( fp0 > fp2 ) { - fMax = fp0; - } else { - fMax = fp2; - } - } else { - if( fp1 > fp2 ) { - fMax = fp1; - } else { - fMax = fp2; - } - } - - // calculate minimum and maximum depth - double fDepthMin = fR - fMin; - double fDepthMax = fMax + fR; - - // if we dont't have overlapping interval - if ( fDepthMin < 0 || fDepthMax < 0 ) { - // do nothing - return false; - } - - double fDepth = 0; - - // if greater depth is on negative side - if ( fDepthMin > fDepthMax ) { - // use smaller depth (one from positive side) - fDepth = fDepthMax; - // flip normal direction - // vNormal[0] = -vNormal[0]; - // vNormal[1] = -vNormal[1]; - // vNormal[2] = -vNormal[2]; - vNormal.scale(-1); - fD = -fD; - // if greater depth is on positive side - } else { - // use smaller depth (one from negative side) - fDepth = fDepthMin; - } - - // if lower depth than best found so far - if (fDepth < m_fBestDepth) { - // remember current axis as best axis - // m_vBestNormal[0] = vNormal[0]; - // m_vBestNormal[1] = vNormal[1]; - // m_vBestNormal[2] = vNormal[2]; - m_vBestNormal.set(vNormal); - m_iBestAxis = iAxis; - //dAASSERT(fDepth>=0); - m_fBestDepth = fDepth; - } - - return true; - } - - // Test cross products of box axis and triangle edges as separating axis - // bool sTrimeshBoxColliderData::_cldTestEdge(dReal fp0, dReal fp1, dReal fR, dReal fD, - // dVector3 vNormal, int iAxis) - boolean _cldTestEdge(double fp0, double fp1, double fR, double fD, - DVector3 vNormal, int iAxis) - { - double fMin, fMax; - - // ===== Begin Patch by Francisco Leon, 2006/10/28 ===== - - // Fixed Null Normal. This prevents boxes passing - // through trimeshes at certain contact angles - - // fMin = vNormal[0] * vNormal[0] + - // vNormal[1] * vNormal[1] + - // vNormal[2] * vNormal[2]; - fMin = vNormal.lengthSquared(); - - if ( fMin <= dEpsilon ) /// THIS NORMAL WOULD BE DANGEROUS - return true; - - // ===== Ending Patch by Francisco Leon ===== - - - // calculate min and max interval values - if ( fp0 < fp1 ) { - fMin = fp0; - fMax = fp1; - } else { - fMin = fp1; - fMax = fp0; - } - - // check if we overlapp - double fDepthMin = fR - fMin; - double fDepthMax = fMax + fR; - - // if we don't overlapp - if ( fDepthMin < 0 || fDepthMax < 0 ) { - // do nothing - return false; - } - - double fDepth; - - // if greater depth is on negative side - if ( fDepthMin > fDepthMax ) { - // use smaller depth (one from positive side) - fDepth = fDepthMax; - // flip normal direction - // vNormal[0] = -vNormal[0]; - // vNormal[1] = -vNormal[1]; - // vNormal[2] = -vNormal[2]; - vNormal.scale(-1); - fD = -fD; - // if greater depth is on positive side - } else { - // use smaller depth (one from negative side) - fDepth = fDepthMin; - } - - // calculate normal's length - double fLength = LENGTHOF(vNormal); - - // if long enough - if ( fLength > 0.0f ) { - - // normalize depth - double fOneOverLength = 1.0f/fLength; - fDepth = fDepth*fOneOverLength; - fD*=fOneOverLength; - - // if lower depth than best found so far (favor face over edges) - if (fDepth*1.5f < m_fBestDepth) { - // remember current axis as best axis - // m_vBestNormal[0] = vNormal[0]*fOneOverLength; - // m_vBestNormal[1] = vNormal[1]*fOneOverLength; - // m_vBestNormal[2] = vNormal[2]*fOneOverLength; - m_vBestNormal.set(vNormal).scale(fOneOverLength); // TODO TZ here positive, above negative, correct? - m_iBestAxis = iAxis; - //dAASSERT(fDepth>=0); - m_fBestDepth = fDepth; - } - } - - return true; - } - - - // clip polygon with plane and generate new polygon points - // static void _cldClipPolyToPlane( dVector3 avArrayIn[], int ctIn, - // dVector3 avArrayOut[], int &ctOut, - // const dVector4 &plPlane ) - void _cldClipPolyToPlane( DVector3 avArrayIn[], int ctIn, - DVector3 avArrayOut[], RefInt ctOut, - final DVector4C plPlane ) - { - // start with no output points - ctOut.i = 0; - - int i0 = ctIn-1; - - // for each edge in input polygon - for (int i1=0; i1= 0 ) { - // emit point - // avArrayOut[ctOut][0] = avArrayIn[i0][0]; - // avArrayOut[ctOut][1] = avArrayIn[i0][1]; - // avArrayOut[ctOut][2] = avArrayIn[i0][2]; - avArrayOut[ctOut.i].set( avArrayIn[i0] ); - ctOut.i++; - } - - // if points are on different sides - if( (fDistance0 > 0 && fDistance1 < 0) || ( fDistance0 < 0 && fDistance1 > 0) ) { - - // find intersection point of edge and plane - DVector3 vIntersectionPoint = new DVector3(); - // vIntersectionPoint[0]= avArrayIn[i0][0] - (avArrayIn[i0][0]-avArrayIn[i1][0])*fDistance0/(fDistance0-fDistance1); - // vIntersectionPoint[1]= avArrayIn[i0][1] - (avArrayIn[i0][1]-avArrayIn[i1][1])*fDistance0/(fDistance0-fDistance1); - // vIntersectionPoint[2]= avArrayIn[i0][2] - (avArrayIn[i0][2]-avArrayIn[i1][2])*fDistance0/(fDistance0-fDistance1); - vIntersectionPoint.eqDiff( avArrayIn[i0], avArrayIn[i1] ); - vIntersectionPoint.scale( - fDistance0/(fDistance0-fDistance1) ); //negate! - vIntersectionPoint.add( avArrayIn[i0] ); - - // emit intersection point - // avArrayOut[ctOut][0] = vIntersectionPoint[0]; - // avArrayOut[ctOut][1] = vIntersectionPoint[1]; - // avArrayOut[ctOut][2] = vIntersectionPoint[2]; - avArrayOut[ctOut.i].set( vIntersectionPoint ); - ctOut.i++; - } - } - - } - - - // bool sTrimeshBoxColliderData::_cldTestSeparatingAxes(const dVector3 &v0, const dVector3 &v1, const dVector3 &v2) { - boolean _cldTestSeparatingAxes(final DVector3C v0, final DVector3C v1, final DVector3C v2) { - // reset best axis - m_iBestAxis = 0; - m_iExitAxis = -1; - m_fBestDepth = MAXVALUE; - - // calculate edges - SUBTRACT(v1,v0,m_vE0); - SUBTRACT(v2,v0,m_vE1); - SUBTRACT(m_vE1,m_vE0,m_vE2); - - // calculate poly normal - //dCROSS(m_vN,=,m_vE0,m_vE1); - m_vN.eqCross(m_vE0, m_vE1); - - // calculate length of face normal - double fNLen = LENGTHOF(m_vN); - - // Even though all triangles might be initially valid, - // a triangle may degenerate into a segment after applying - // space transformation. - if (fNLen == 0) { - return false; - } - - // extract box axes as vectors - DVector3 vA0 = new DVector3(), vA1 = new DVector3(), vA2 = new DVector3(); - GETCOL(m_mHullBoxRot,0,vA0); - GETCOL(m_mHullBoxRot,1,vA1); - GETCOL(m_mHullBoxRot,2,vA2); - - // box halfsizes - double fa0 = m_vBoxHalfSize.get0(); - double fa1 = m_vBoxHalfSize.get1(); - double fa2 = m_vBoxHalfSize.get2(); - - // calculate relative position between box and triangle - DVector3 vD = new DVector3(); - SUBTRACT(v0,m_vHullBoxPos,vD); - - DVector3 vL = new DVector3(); - double fp0, fp1, fp2, fR, fD; - - // Test separating axes for intersection - // ************************************************ - // Axis 1 - Triangle Normal - SET(vL,m_vN); - fp0 = vL.dot(vD); - fp1 = fp0; - fp2 = fp0; - fR=fa0*dFabs( m_vN.dot(vA0) ) + fa1 * dFabs( m_vN.dot(vA1) ) + fa2 * dFabs( m_vN.dot(vA2) ); - - if (!_cldTestNormal(fp0, fR, vL, 1)) { - m_iExitAxis=1; - return false; - } - - // ************************************************ - - // Test Faces - // ************************************************ - // Axis 2 - Box X-Axis - SET(vL,vA0); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 + vA0.dot(m_vE0); - fp2 = fp0 + vA0.dot(m_vE1); - fR = fa0; - - if (!_cldTestFace(fp0, fp1, fp2, fR, fD, vL, 2)) { - m_iExitAxis=2; - return false; - } - // ************************************************ - - // ************************************************ - // Axis 3 - Box Y-Axis - SET(vL,vA1); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 + vA1.dot(m_vE0); - fp2 = fp0 + vA1.dot(m_vE1); - fR = fa1; - - if (!_cldTestFace(fp0, fp1, fp2, fR, fD, vL, 3)) { - m_iExitAxis=3; - return false; - } - - // ************************************************ - - // ************************************************ - // Axis 4 - Box Z-Axis - SET(vL,vA2); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 + vA2.dot(m_vE0); - fp2 = fp0 + vA2.dot(m_vE1); - fR = fa2; - - if (!_cldTestFace(fp0, fp1, fp2, fR, fD, vL, 4)) { - m_iExitAxis=4; - return false; - } - - // ************************************************ - - // Test Edges - // ************************************************ - // Axis 5 - Box X-Axis cross Edge0 - //dCROSS(vL,=,vA0,m_vE0); - vL.eqCross(vA0, m_vE0); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0; - fp2 = fp0 + vA0.dot(m_vN); - fR = fa1 * dFabs(vA2.dot(m_vE0)) + fa2 * dFabs(vA1.dot(m_vE0)); - - if (!_cldTestEdge(fp1, fp2, fR, fD, vL, 5)) { - m_iExitAxis=5; - return false; - } - // ************************************************ - - // ************************************************ - // Axis 6 - Box X-Axis cross Edge1 - //dCROSS(vL,=,vA0,m_vE1); - vL.eqCross(vA0, m_vE1); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 - vA0.dot(m_vN); - fp2 = fp0; - fR = fa1 * dFabs(vA2.dot(m_vE1)) + fa2 * dFabs(vA1.dot(m_vE1)); - - if (!_cldTestEdge(fp0, fp1, fR, fD, vL, 6)) { - m_iExitAxis=6; - return false; - } - // ************************************************ - - // ************************************************ - // Axis 7 - Box X-Axis cross Edge2 - //dCROSS(vL,=,vA0,m_vE2); - vL.eqCross(vA0, m_vE2); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 - vA0.dot(m_vN); - fp2 = fp0 - vA0.dot(m_vN); - fR = fa1 * dFabs(vA2.dot(m_vE2)) + fa2 * dFabs(vA1.dot(m_vE2)); - - if (!_cldTestEdge(fp0, fp1, fR, fD, vL, 7)) { - m_iExitAxis=7; - return false; - } - - // ************************************************ - - // ************************************************ - // Axis 8 - Box Y-Axis cross Edge0 - //dCROSS(vL,=,vA1,m_vE0); - vL.eqCross(vA1, m_vE0); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0; - fp2 = fp0 + vA1.dot(m_vN); - fR = fa0 * dFabs(vA2.dot(m_vE0)) + fa2 * dFabs(vA0.dot(m_vE0)); - - if (!_cldTestEdge(fp0, fp2, fR, fD, vL, 8)) { - m_iExitAxis=8; - return false; - } - - // ************************************************ - - // ************************************************ - // Axis 9 - Box Y-Axis cross Edge1 - //dCROSS(vL,=,vA1,m_vE1); - vL.eqCross(vA1, m_vE1); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 - vA1.dot(m_vN); - fp2 = fp0; - fR = fa0 * dFabs(vA2.dot(m_vE1)) + fa2 * dFabs(vA0.dot(m_vE1)); - - if (!_cldTestEdge(fp0, fp1, fR, fD, vL, 9)) { - m_iExitAxis=9; - return false; - } - - // ************************************************ - - // ************************************************ - // Axis 10 - Box Y-Axis cross Edge2 - //dCROSS(vL,=,vA1,m_vE2); - vL.eqCross(vA1, m_vE2); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 - vA1.dot(m_vN); - fp2 = fp0 - vA1.dot(m_vN); - fR = fa0 * dFabs(vA2.dot(m_vE2)) + fa2 * dFabs(vA0.dot(m_vE2)); - - if (!_cldTestEdge(fp0, fp1, fR, fD, vL, 10)) { - m_iExitAxis=10; - return false; - } - - // ************************************************ - - // ************************************************ - // Axis 11 - Box Z-Axis cross Edge0 - //dCROSS(vL,=,vA2,m_vE0); - vL.eqCross(vA2, m_vE0); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0; - fp2 = fp0 + vA2.dot(m_vN); - fR = fa0 * dFabs(vA1.dot(m_vE0)) + fa1 * dFabs(vA0.dot(m_vE0)); - - if (!_cldTestEdge(fp0, fp2, fR, fD, vL, 11)) { - m_iExitAxis=11; - return false; - } - // ************************************************ - - // ************************************************ - // Axis 12 - Box Z-Axis cross Edge1 - //dCROSS(vL,=,vA2,m_vE1); - vL.eqCross(vA2, m_vE1); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 - vA2.dot(m_vN); - fp2 = fp0; - fR = fa0 * dFabs(vA1.dot(m_vE1)) + fa1 * dFabs(vA0.dot(m_vE1)); - - if (!_cldTestEdge(fp0, fp1, fR, fD, vL, 12)) { - m_iExitAxis=12; - return false; - } - // ************************************************ - - // ************************************************ - // Axis 13 - Box Z-Axis cross Edge2 - // OdeMath.dCROSS(vL,OP.EQ,vA2,m_vE2); - vL.eqCross(vA2, m_vE2); - fD = vL.dot(m_vN)/fNLen; - fp0 = vL.dot(vD); - fp1 = fp0 - vA2.dot(m_vN); - fp2 = fp0 - vA2.dot(m_vN); - fR = fa0 * dFabs(vA1.dot(m_vE2)) + fa1 * dFabs(vA0.dot(m_vE2)); - - if (!_cldTestEdge(fp0, fp1, fR, fD, vL, 13)) { - m_iExitAxis=13; - return false; - } - - // ************************************************ - return true; - } - - - - - - // find two closest points on two lines - // static bool _cldClosestPointOnTwoLines( dVector3 vPoint1, dVector3 vLenVec1, - // dVector3 vPoint2, dVector3 vLenVec2, - // dReal &fvalue1, dReal &fvalue2) - boolean _cldClosestPointOnTwoLines( DVector3C vPoint1, DVector3C vLenVec1, - DVector3C vPoint2, DVector3C vLenVec2, - RefDouble fvalue1, RefDouble fvalue2) - { - // calculate denominator - DVector3 vp = new DVector3(); - SUBTRACT(vPoint2,vPoint1,vp); - double fuaub = vLenVec1.dot(vLenVec2); - double fq1 = vLenVec1.dot(vp); - double fq2 = -vLenVec2.dot(vp); - double fd = 1.0f - fuaub * fuaub; - - // if denominator is positive - if (fd > 0.0f) { - // calculate points of closest approach - fd = 1.0f/fd; - fvalue1.d = (fq1 + fuaub*fq2)*fd; - fvalue2.d = (fuaub*fq1 + fq2)*fd; - return true; - // otherwise - } else { - // lines are parallel - fvalue1.d = 0.0f; - fvalue2.d = 0.0f; - return false; - } - } - - - - - - // clip and generate contacts - // void sTrimeshBoxColliderData::_cldClipping(const dVector3 &v0, const dVector3 &v1, const dVector3 &v2, int TriIndex) { - void _cldClipping(final DVector3C v0, final DVector3C v1, final DVector3C v2, final int TriIndex) { - dIASSERT( ((m_iFlags & CONTACTS_UNIMPORTANT)==0) || m_ctContacts.i < (m_iFlags & DxGeom.NUMC_MASK) ); // Do not call the function if there is no room to store results - - // if we have edge/edge intersection - if (m_iBestAxis > 4 ) { - DVector3 vub = new DVector3(), vPb = new DVector3(), vPa = new DVector3(); - - SET(vPa,m_vHullBoxPos); - - // calculate point on box edge - for( int i=0; i<3; i++) { - DVector3 vRotCol = new DVector3(); - GETCOL(m_mHullBoxRot,i,vRotCol); - double fSign = m_vBestNormal.dot(vRotCol) > 0 ? 1.0f : -1.0f; - - // vPa[0] += fSign * m_vBoxHalfSize[i] * vRotCol[0]; - // vPa[1] += fSign * m_vBoxHalfSize[i] * vRotCol[1]; - // vPa[2] += fSign * m_vBoxHalfSize[i] * vRotCol[2]; - vPa.addScaled( vRotCol, fSign * m_vBoxHalfSize.get(i)); //TODO use colView! - } - - int iEdge = (m_iBestAxis-5)%3; - - // decide which edge is on triangle - if ( iEdge == 0 ) { - SET(vPb,v0); - SET(vub,m_vE0); - } else if ( iEdge == 1) { - SET(vPb,v2); - SET(vub,m_vE1); - } else { - SET(vPb,v1); - SET(vub,m_vE2); - } - - - // setup direction parameter for face edge - //dNormalize3(vub); - vub.normalize(); - - RefDouble fParam1 = new RefDouble(), fParam2 = new RefDouble(); - - // setup direction parameter for box edge - DVector3 vua = new DVector3(); - int col=(m_iBestAxis-5)/3; - GETCOL(m_mHullBoxRot,col,vua); - - // find two closest points on both edges - _cldClosestPointOnTwoLines( vPa, vua, vPb, vub, fParam1, fParam2 ); - // vPa[0] += vua[0]*fParam1; - // vPa[1] += vua[1]*fParam1; - // vPa[2] += vua[2]*fParam1; - vPa.addScaled( vua, fParam1.d ); - - // vPb[0] += vub[0]*fParam2; - // vPb[1] += vub[1]*fParam2; - // vPb[2] += vub[2]*fParam2; - vPb.addScaled( vub, fParam2.d ); - - // calculate collision point - DVector3 vPntTmp = new DVector3(); - ADD(vPa,vPb,vPntTmp); - - // vPntTmp[0]*=0.5f; - // vPntTmp[1]*=0.5f; - // vPntTmp[2]*=0.5f; - vPntTmp.scale( 0.5 ); - - // generate contact point between two closest points - // #if 0 //#ifdef ORIG -- if to use conditional define, GenerateContact must be moved into #else - // dContactGeom* Contact = SAFECONTACT(m_iFlags, m_ContactGeoms, m_ctContacts, m_iStride); - // Contact->depth = m_fBestDepth; - // SET(Contact->normal,m_vBestNormal); - // SET(Contact->pos,vPntTmp); - // Contact->g1 = Geom1; - // Contact->g2 = Geom2; - // Contact->side1 = TriIndex; - // Contact->side2 = -1; - // m_ctContacts++; - // #endif - GenerateContact(m_iFlags, m_ContactGeoms, m_iStride, m_Geom1, m_Geom2, TriIndex, - vPntTmp, m_vBestNormal, m_fBestDepth, m_ctContacts); - - - // if triangle is the referent face then clip box to triangle face - } else if (m_iBestAxis == 1) { - - DVector3 vNormal2 = new DVector3(); - // vNormal2[0]=-m_vBestNormal[0]; - // vNormal2[1]=-m_vBestNormal[1]; - // vNormal2[2]=-m_vBestNormal[2]; - vNormal2.set( m_vBestNormal ).scale( -1 ); - - - // vNr is normal in box frame, pointing from triangle to box - DMatrix3 mTransposed;// = new DMatrix3(); - // mTransposed[0*4+0]=m_mHullBoxRot[0*4+0]; - // mTransposed[0*4+1]=m_mHullBoxRot[1*4+0]; - // mTransposed[0*4+2]=m_mHullBoxRot[2*4+0]; - // - // mTransposed[1*4+0]=m_mHullBoxRot[0*4+1]; - // mTransposed[1*4+1]=m_mHullBoxRot[1*4+1]; - // mTransposed[1*4+2]=m_mHullBoxRot[2*4+1]; - // - // mTransposed[2*4+0]=m_mHullBoxRot[0*4+2]; - // mTransposed[2*4+1]=m_mHullBoxRot[1*4+2]; - // mTransposed[2*4+2]=m_mHullBoxRot[2*4+2]; - mTransposed = m_mHullBoxRot.reTranspose(); - - DVector3 vNr = new DVector3(); - // vNr[0]=mTransposed[0*4+0]*vNormal2[0]+ mTransposed[0*4+1]*vNormal2[1]+ mTransposed[0*4+2]*vNormal2[2]; - // vNr[1]=mTransposed[1*4+0]*vNormal2[0]+ mTransposed[1*4+1]*vNormal2[1]+ mTransposed[1*4+2]*vNormal2[2]; - // vNr[2]=mTransposed[2*4+0]*vNormal2[0]+ mTransposed[2*4+1]*vNormal2[1]+ mTransposed[2*4+2]*vNormal2[2]; - vNr.eqProd(mTransposed, vNormal2); - - - DVector3 vAbsNormal = new DVector3(); - // vAbsNormal[0] = dFabs( vNr[0] ); - // vAbsNormal[1] = dFabs( vNr[1] ); - // vAbsNormal[2] = dFabs( vNr[2] ); - vAbsNormal.set( vNr ).eqAbs(); - - // get closest face from box - int iB0, iB1, iB2; - if (vAbsNormal.get1() > vAbsNormal.get0()) { - if (vAbsNormal.get1() > vAbsNormal.get2()) { - iB1 = 0; iB0 = 1; iB2 = 2; - } else { - iB1 = 0; iB2 = 1; iB0 = 2; - } - } else { - - if (vAbsNormal.get0() > vAbsNormal.get2()) { - iB0 = 0; iB1 = 1; iB2 = 2; - } else { - iB1 = 0; iB2 = 1; iB0 = 2; - } - } - - // Here find center of box face we are going to project - DVector3 vCenter = new DVector3(); - DVector3 vRotCol = new DVector3(); - GETCOL(m_mHullBoxRot,iB0,vRotCol); - - if (vNr.get(iB0) > 0) { - // vCenter[0] = m_vHullBoxPos[0] - v0[0] - m_vBoxHalfSize[iB0] * vRotCol[0]; - // vCenter[1] = m_vHullBoxPos[1] - v0[1] - m_vBoxHalfSize[iB0] * vRotCol[1]; - // vCenter[2] = m_vHullBoxPos[2] - v0[2] - m_vBoxHalfSize[iB0] * vRotCol[2]; - vCenter.eqSum(v0, -1, vRotCol, -m_vBoxHalfSize.get(iB0)); - vCenter.add(m_vHullBoxPos); - } else { - // vCenter[0] = m_vHullBoxPos[0] - v0[0] + m_vBoxHalfSize[iB0] * vRotCol[0]; - // vCenter[1] = m_vHullBoxPos[1] - v0[1] + m_vBoxHalfSize[iB0] * vRotCol[1]; - // vCenter[2] = m_vHullBoxPos[2] - v0[2] + m_vBoxHalfSize[iB0] * vRotCol[2]; - vCenter.eqSum(v0, -1, vRotCol, m_vBoxHalfSize.get(iB0)); - vCenter.add(m_vHullBoxPos); - } - - // Here find 4 corner points of box - DVector3[] avPoints = { new DVector3(), new DVector3(), new DVector3(), new DVector3() };//[4]; - - DVector3 vRotCol2 = new DVector3(); - GETCOL(m_mHullBoxRot,iB1,vRotCol); - GETCOL(m_mHullBoxRot,iB2,vRotCol2); - - // for(int x=0;x<3;x++) { - // avPoints[0][x] = vCenter[x] + (m_vBoxHalfSize[iB1] * vRotCol[x]) - (m_vBoxHalfSize[iB2] * vRotCol2[x]); - // avPoints[1][x] = vCenter[x] - (m_vBoxHalfSize[iB1] * vRotCol[x]) - (m_vBoxHalfSize[iB2] * vRotCol2[x]); - // avPoints[2][x] = vCenter[x] - (m_vBoxHalfSize[iB1] * vRotCol[x]) + (m_vBoxHalfSize[iB2] * vRotCol2[x]); - // avPoints[3][x] = vCenter[x] + (m_vBoxHalfSize[iB1] * vRotCol[x]) + (m_vBoxHalfSize[iB2] * vRotCol2[x]); - // } - double tz1 = m_vBoxHalfSize.get(iB1); - double tz2 = m_vBoxHalfSize.get(iB2); - avPoints[0].eqSum( vRotCol, tz1, vRotCol2, -tz2).add (vCenter); - avPoints[1].eqSum( vRotCol, -tz1, vRotCol2, -tz2).add (vCenter); - avPoints[2].eqSum( vRotCol, -tz1, vRotCol2, tz2).add (vCenter); - avPoints[3].eqSum( vRotCol, tz1, vRotCol2, tz2).add (vCenter); - - // clip Box face with 4 planes of triangle (1 face plane, 3 egde planes) - DVector3[] avTempArray1 = DVector3.newArray(9); - DVector3[] avTempArray2 = DVector3.newArray(9); - DVector4 plPlane = new DVector4(); - - RefInt iTempCnt1 = new RefInt(0); - RefInt iTempCnt2 = new RefInt(0); - - // zeroify vectors - necessary? - // for(int i=0; i<9; i++) { - // // avTempArray1[i][0]=0; - // // avTempArray1[i][1]=0; - // // avTempArray1[i][2]=0; - // avTempArray1[i] = new DVector3(); - // - // // avTempArray2[i][0]=0; - // // avTempArray2[i][1]=0; - // // avTempArray2[i][2]=0; - // avTempArray2[i] = new DVector3(); - // } - - - // Normal plane - DVector3 vTemp = new DVector3(); - // vTemp[0]=-m_vN[0]; - // vTemp[1]=-m_vN[1]; - // vTemp[2]=-m_vN[2]; - vTemp.set( m_vN ).scale(-1); - //dNormalize3(vTemp); - vTemp.normalize(); - CONSTRUCTPLANE(plPlane,vTemp,0); - - _cldClipPolyToPlane( avPoints, 4, avTempArray1, iTempCnt1, plPlane ); - - - // Plane p0 - DVector3 vTemp2 = new DVector3(); - SUBTRACT(v1,v0,vTemp2); - //dCROSS(vTemp,=,m_vN,vTemp2); - vTemp.eqCross(m_vN, vTemp2); - //dNormalize3(vTemp); - vTemp.normalize(); - CONSTRUCTPLANE(plPlane,vTemp,0); - - _cldClipPolyToPlane( avTempArray1, iTempCnt1.i, avTempArray2, iTempCnt2, plPlane ); - - // Plane p1 - SUBTRACT(v2,v1,vTemp2); - //dCROSS(vTemp,=,m_vN,vTemp2); - vTemp.eqCross(m_vN, vTemp2); - //dNormalize3(vTemp); - vTemp.normalize(); - SUBTRACT(v0,v2,vTemp2); - CONSTRUCTPLANE(plPlane,vTemp,vTemp2.dot(vTemp)); - - _cldClipPolyToPlane( avTempArray2, iTempCnt2.i, avTempArray1, iTempCnt1, plPlane ); - - // Plane p2 - SUBTRACT(v0,v2,vTemp2); - //dCROSS(vTemp,=,m_vN,vTemp2); - vTemp.eqCross(m_vN, vTemp2); - //dNormalize3(vTemp); - vTemp.normalize(); - CONSTRUCTPLANE(plPlane,vTemp,0); - - _cldClipPolyToPlane( avTempArray1, iTempCnt1.i, avTempArray2, iTempCnt2, plPlane ); - - // END of clipping polygons - - // for each generated contact point - for ( int i=0; i 0) { - fTempDepth = 0; - } - - DVector3 vPntTmp = new DVector3(); - ADD(avTempArray2[i],v0,vPntTmp); - - // #if 0 //#ifdef ORIG -- if to use conditional define, GenerateContact must be moved into #else - // dContactGeom* Contact = SAFECONTACT(m_iFlags, m_ContactGeoms, m_ctContacts, m_iStride); - // - // Contact->depth = -fTempDepth; - // SET(Contact->normal,m_vBestNormal); - // SET(Contact->pos,vPntTmp); - // Contact->g1 = Geom1; - // Contact->g2 = Geom2; - // Contact->side1 = TriIndex; - // Contact->side2 = -1; - // m_ctContacts++; - // #endif - GenerateContact(m_iFlags, m_ContactGeoms, m_iStride, m_Geom1, m_Geom2, TriIndex, - vPntTmp, m_vBestNormal, -fTempDepth, m_ctContacts); - - if ((m_ctContacts.i | CONTACTS_UNIMPORTANT) == (m_iFlags & (DxGeom.NUMC_MASK | CONTACTS_UNIMPORTANT))) { - break; - } - } - - //dAASSERT(m_ctContacts>0); - - // if box face is the referent face, then clip triangle on box face - } else { // 2 <= if iBestAxis <= 4 - - // get normal of box face - DVector3 vNormal2 = new DVector3(); - SET(vNormal2,m_vBestNormal); - - // get indices of box axes in correct order - int iA0,iA1,iA2; - iA0 = m_iBestAxis-2; - if ( iA0 == 0 ) { - iA1 = 1; iA2 = 2; - } else if ( iA0 == 1 ) { - iA1 = 0; iA2 = 2; - } else { - iA1 = 0; iA2 = 1; - } - - DVector3[] avPoints = { new DVector3(), new DVector3(), new DVector3() };//new DVector3[3]; - // calculate triangle vertices in box frame - SUBTRACT(v0,m_vHullBoxPos,avPoints[0]); - SUBTRACT(v1,m_vHullBoxPos,avPoints[1]); - SUBTRACT(v2,m_vHullBoxPos,avPoints[2]); - - // CLIP Polygons - // define temp data for clipping - DVector3[] avTempArray1 = DVector3.newArray(9); - DVector3[] avTempArray2 = DVector3.newArray(9); - - RefInt iTempCnt1 = new RefInt(), iTempCnt2 = new RefInt(); - - // zeroify vectors - necessary? - // for(int i=0; i<9; i++) { - // // avTempArray1[i][0]=0; - // // avTempArray1[i][1]=0; - // // avTempArray1[i][2]=0; - // avTempArray1[i] = new DVector3(); - // - // // avTempArray2[i][0]=0; - // // avTempArray2[i][1]=0; - // // avTempArray2[i][2]=0; - // avTempArray2[i] = new DVector3(); - // } - - // clip triangle with 5 box planes (1 face plane, 4 edge planes) - - DVector4 plPlane = new DVector4(); - - // Normal plane - DVector3 vTemp = new DVector3(); - // vTemp[0]=-vNormal2[0]; - // vTemp[1]=-vNormal2[1]; - // vTemp[2]=-vNormal2[2]; - vTemp.set(vNormal2).scale(-1); - CONSTRUCTPLANE(plPlane,vTemp,m_vBoxHalfSize.get(iA0)); - - _cldClipPolyToPlane( avPoints, 3, avTempArray1, iTempCnt1, plPlane ); - - - // Plane p0 - GETCOL(m_mHullBoxRot,iA1,vTemp); - CONSTRUCTPLANE(plPlane,vTemp,m_vBoxHalfSize.get(iA1)); - - _cldClipPolyToPlane( avTempArray1, iTempCnt1.i, avTempArray2, iTempCnt2, plPlane ); - - - // Plane p1 - GETCOL(m_mHullBoxRot,iA1,vTemp); - // vTemp[0]=-vTemp[0]; - // vTemp[1]=-vTemp[1]; - // vTemp[2]=-vTemp[2]; - vTemp.scale(-1); - CONSTRUCTPLANE(plPlane,vTemp,m_vBoxHalfSize.get(iA1)); - - _cldClipPolyToPlane( avTempArray2, iTempCnt2.i, avTempArray1, iTempCnt1, plPlane ); - - // Plane p2 - GETCOL(m_mHullBoxRot,iA2,vTemp); - CONSTRUCTPLANE(plPlane,vTemp,m_vBoxHalfSize.get(iA2)); - - _cldClipPolyToPlane( avTempArray1, iTempCnt1.i, avTempArray2, iTempCnt2, plPlane ); - - // Plane p3 - GETCOL(m_mHullBoxRot,iA2,vTemp); - // vTemp[0]=-vTemp[0]; - // vTemp[1]=-vTemp[1]; - // vTemp[2]=-vTemp[2]; - vTemp.scale(-1); - CONSTRUCTPLANE(plPlane,vTemp,m_vBoxHalfSize.get(iA2)); - - _cldClipPolyToPlane( avTempArray2, iTempCnt2.i, avTempArray1, iTempCnt1, plPlane ); - - - // for each generated contact point - for ( int i=0; i 0) { - fTempDepth = 0; - } - - // generate contact data - DVector3 vPntTmp = new DVector3(); - ADD(avTempArray1[i],m_vHullBoxPos,vPntTmp); - - // #if 0 //#ifdef ORIG -- if to use conditional define, GenerateContact must be moved into #else - // dContactGeom* Contact = SAFECONTACT(m_iFlags, m_ContactGeoms, m_ctContacts, m_iStride); - // - // Contact->depth = -fTempDepth; - // SET(Contact->normal,m_vBestNormal); - // SET(Contact->pos,vPntTmp); - // Contact->g1 = Geom1; - // Contact->g2 = Geom2; - // Contact->side1 = TriIndex; - // Contact->side2 = -1; - // m_ctContacts++; - // #endif - GenerateContact(m_iFlags, m_ContactGeoms, m_iStride, m_Geom1, m_Geom2, TriIndex, - vPntTmp, m_vBestNormal, -fTempDepth, m_ctContacts); - - if ((m_ctContacts.i | CONTACTS_UNIMPORTANT) == (m_iFlags & (DxGeom.NUMC_MASK | CONTACTS_UNIMPORTANT))) { - break; - } - } - - //dAASSERT(m_ctContacts>0); - } - } - - - - - - // test one mesh triangle on intersection with given box - // void sTrimeshBoxColliderData::_cldTestOneTriangle(const dVector3 &v0, const dVector3 &v1, const dVector3 &v2, int TriIndex)//, void *pvUser) - private void _cldTestOneTriangle(final DVector3C v0, final DVector3C v1, - final DVector3C v2, int TriIndex)//, void *pvUser) - { - // do intersection test and find best separating axis - if(!_cldTestSeparatingAxes(v0, v1, v2)) { - // if not found do nothing - return; - } - - // if best separation axis is not found - if (m_iBestAxis == 0) { - // this should not happen (we should already exit in that case) - //dMessage (0, "best separation axis not found"); - // do nothing - return; - } - - _cldClipping(v0, v1, v2, TriIndex); - } - - - // void sTrimeshBoxColliderData::SetupInitialContext(dxTriMesh *TriMesh, dxGeom *BoxGeom, - // int Flags, dContactGeom* Contacts, int Stride) - void SetupInitialContext(DxTriMesh TriMesh, DxBox BoxGeom, - int Flags, DContactGeomBuffer Contacts, int Stride) - { - // // get source hull position, orientation and half size - // const dMatrix3& mRotBox=*(const dMatrix3*)dGeomGetRotation(BoxGeom); - // const dVector3& vPosBox=*(const dVector3*)dGeomGetPosition(BoxGeom); - // - // // to global - // SETM(m_mHullBoxRot,mRotBox); - // SET(m_vHullBoxPos,vPosBox); - m_mHullBoxRot.set(BoxGeom.getRotation()); - m_vHullBoxPos.set(BoxGeom.getPosition()); - - //dGeomBoxGetLengths(BoxGeom, m_vBoxHalfSize); - m_vBoxHalfSize.set(BoxGeom.getLengths()); - // m_vBoxHalfSize[0] *= 0.5f; - // m_vBoxHalfSize[1] *= 0.5f; - // m_vBoxHalfSize[2] *= 0.5f; - m_vBoxHalfSize.scale(0.5f); - - // get destination hull position and orientation - // const dVector3& vPosMesh=*(const dVector3*)dGeomGetPosition(TriMesh); - // - // // to global - // SET(m_vHullDstPos,vPosMesh); - m_vHullDstPos.set(TriMesh.getPosition()); - - // global info for contact creation - m_ctContacts.i = 0; - if (Stride != 1) throw new IllegalArgumentException("stride = " + Stride); - m_iStride=Stride; - m_iFlags=Flags; - m_ContactGeoms=Contacts; - m_Geom1=TriMesh; - m_Geom2=BoxGeom; - - // reset stuff - m_fBestDepth = MAXVALUE; - // m_vBestNormal[0]=0; - // m_vBestNormal[1]=0; - // m_vBestNormal[2]=0; - m_vBestNormal.setZero(); - } - - // int sTrimeshBoxColliderData::TestCollisionForSingleTriangle(int ctContacts0, int Triint, - // dVector3 dv[3], bool &bOutFinishSearching) - private int TestCollisionForSingleTriangle(int ctContacts0, int Triint, - DVector3 dv[], RefBoolean bOutFinishSearching) - { - // test this triangle - _cldTestOneTriangle(dv[0],dv[1],dv[2],Triint); - - // fill-in tri index for generated contacts - for (; ctContacts0 < m_ctContacts.i; ctContacts0++) { - //DContactGeom pContact = SAFECONTACT(m_iFlags, m_ContactGeoms, ctContacts0, m_iStride); - DContactGeom pContact = m_ContactGeoms.getSafe(m_iFlags, ctContacts0); - pContact.side1 = Triint; - pContact.side2 = -1; - } - - /* - NOTE by Oleh_Derevenko: - The function continues checking triangles after maximal number - of contacts is reached because it selects maximal penetration depths. - See also comments in GenerateContact() - */ - bOutFinishSearching.b = ((m_ctContacts.i | CONTACTS_UNIMPORTANT) == (m_iFlags & (DxGeom.NUMC_MASK | CONTACTS_UNIMPORTANT))); - - return ctContacts0; - } - - - // // OPCODE version of box to mesh collider - // #if dTRIMESH_OPCODE - // static void dQueryBTLPotentialCollisionTriangles(OBBCollider &Collider, - // const sTrimeshBoxColliderData &cData, dxTriMesh *TriMesh, dxGeom *BoxGeom, - // OBBCache &BoxCache) - // { - // // get source hull position, orientation and half size - // const dMatrix3& mRotBox=*(const dMatrix3*)dGeomGetRotation(BoxGeom); - // const dVector3& vPosBox=*(const dVector3*)dGeomGetPosition(BoxGeom); - // - // // Make OBB - // OBB Box; - // Box.mCenter.x = vPosBox[0]; - // Box.mCenter.y = vPosBox[1]; - // Box.mCenter.z = vPosBox[2]; - // - // // It is a potential issue to explicitly cast to float - // // if custom width floating point type is introduced in OPCODE. - // // It is necessary to make a typedef and cast to it - // // (e.g. typedef float opc_float;) - // // However I'm not sure in what header it should be added. - // - // Box.mExtents.x = /*(float)*/cData.m_vBoxHalfSize[0]; - // Box.mExtents.y = /*(float)*/cData.m_vBoxHalfSize[1]; - // Box.mExtents.z = /*(float)*/cData.m_vBoxHalfSize[2]; - // - // Box.mRot.m[0][0] = /*(float)*/mRotBox[0]; - // Box.mRot.m[1][0] = /*(float)*/mRotBox[1]; - // Box.mRot.m[2][0] = /*(float)*/mRotBox[2]; - // - // Box.mRot.m[0][1] = /*(float)*/mRotBox[4]; - // Box.mRot.m[1][1] = /*(float)*/mRotBox[5]; - // Box.mRot.m[2][1] = /*(float)*/mRotBox[6]; - // - // Box.mRot.m[0][2] = /*(float)*/mRotBox[8]; - // Box.mRot.m[1][2] = /*(float)*/mRotBox[9]; - // Box.mRot.m[2][2] = /*(float)*/mRotBox[10]; - // - // Matrix4x4 amatrix; - // Matrix4x4 BoxMatrix = MakeMatrix(vPosBox, mRotBox, amatrix); - // - // Matrix4x4 InvBoxMatrix; - // InvertPRMatrix(InvBoxMatrix, BoxMatrix); - // - // // get destination hull position and orientation - // const dMatrix3& mRotMesh=*(const dMatrix3*)dGeomGetRotation(TriMesh); - // const dVector3& vPosMesh=*(const dVector3*)dGeomGetPosition(TriMesh); - // - // // TC results - // if (TriMesh->doBoxTC) { - // dxTriMesh::BoxTC* BoxTC = 0; - // for (int i = 0; i < TriMesh->BoxTCCache.size(); i++){ - // if (TriMesh->BoxTCCache[i].Geom == BoxGeom){ - // BoxTC = &TriMesh->BoxTCCache[i]; - // break; - // } - // } - // if (!BoxTC){ - // TriMesh->BoxTCCache.push(dxTriMesh::BoxTC()); - // - // BoxTC = &TriMesh->BoxTCCache[TriMesh->BoxTCCache.size() - 1]; - // BoxTC->Geom = BoxGeom; - // BoxTC->FatCoeff = 1.1f; // Pierre recommends this, instead of 1.0 - // } - // - // // Intersect - // Collider.SetTemporalCoherence(true); - // Collider.Collide(*BoxTC, Box, TriMesh->Data->BVTree, null, &MakeMatrix(vPosMesh, mRotMesh, amatrix)); - // } - // else { - // Collider.SetTemporalCoherence(false); - // Collider.Collide(BoxCache, Box, TriMesh->Data->BVTree, null, - // &MakeMatrix(vPosMesh, mRotMesh, amatrix)); - // } - // } - // - // int dCollideBTL(dxGeom* g1, dxGeom* BoxGeom, int Flags, dContactGeom* Contacts, int Stride){ - // dIASSERT (Stride >= (int)sizeof(dContactGeom)); - // dIASSERT (g1->type == dTriMeshClass); - // dIASSERT (BoxGeom->type == dBoxClass); - // dIASSERT ((Flags & NUMC_MASK) >= 1); - // - // dxTriMesh* TriMesh = (dxTriMesh*)g1; - // - // sTrimeshBoxColliderData cData; - // cData.SetupInitialContext(TriMesh, BoxGeom, Flags, Contacts, Stride); - // - // const unsigned uiTLSKind = TriMesh->getParentSpaceTLSKind(); - // dIASSERT(uiTLSKind == BoxGeom->getParentSpaceTLSKind()); // The colliding spaces must use matching cleanup method - // TrimeshCollidersCache *pccColliderCache = GetTrimeshCollidersCache(uiTLSKind); - // OBBCollider& Collider = pccColliderCache->_OBBCollider; - // - // dQueryBTLPotentialCollisionTriangles(Collider, cData, TriMesh, BoxGeom, - // pccColliderCache->defaultBoxCache); - // - // if (!Collider.GetContactStatus()) { - // // no collision occurred - // return 0; - // } - // - // // Retrieve data - // int TriCount = Collider.GetNbTouchedPrimitives(); - // const int* Triangles = (const int*)Collider.GetTouchedPrimitives(); - // - // if (TriCount != 0){ - // if (TriMesh->ArrayCallback != null){ - // TriMesh->ArrayCallback(TriMesh, BoxGeom, Triangles, TriCount); - // } - // - // // get destination hull position and orientation - // const dMatrix3& mRotMesh=*(const dMatrix3*)dGeomGetRotation(TriMesh); - // const dVector3& vPosMesh=*(const dVector3*)dGeomGetPosition(TriMesh); - // - // int ctContacts0 = 0; - // - // // loop through all intersecting triangles - // for (int i = 0; i < TriCount; i++){ - // const int Triint = Triangles[i]; - // if (!Callback(TriMesh, BoxGeom, Triint)) continue; - // - // dVector3 dv[3]; - // FetchTriangle(TriMesh, Triint, vPosMesh, mRotMesh, dv); - // - // bool bFinishSearching; - // ctContacts0 = cData.TestCollisionForSingleTriangle(ctContacts0, Triint, dv, bFinishSearching); - // - // if (bFinishSearching) { - // break; - // } - // } - // } - // - // return cData.m_ctContacts; - // } - // #endif - - } //TZ end of data? - - // GIMPACT version of box to mesh collider - // #if dTRIMESH_GIMPACT - // int dCollideBTL(dxGeom* g1, dxGeom* BoxGeom, int Flags, dContactGeom* Contacts, int Stride) - int dCollideBTL(DxTriMesh g1, DxBox BoxGeom, int Flags, DContactGeomBuffer Contacts, int Stride) - { - dIASSERT (Stride >= 1);//(int)sizeof(dContactGeom)); - // dIASSERT (g1.type == dTriMeshClass); - // dIASSERT (BoxGeom.type == dBoxClass); - dIASSERT ((Flags & DxGeom.NUMC_MASK) >= 1); - - - DxGimpact TriMesh = (DxGimpact) g1; - - g1.recomputeAABB(); - BoxGeom.recomputeAABB(); - - - sTrimeshBoxColliderData cData = new sTrimeshBoxColliderData(); - cData.SetupInitialContext(TriMesh, BoxGeom, Flags, Contacts, Stride); - - //*****at first , collide box aabb******// - - //GIM_TRIMESH * ptrimesh = &TriMesh.m_collision_trimesh; - GimTrimesh ptrimesh = TriMesh.m_collision_trimesh(); - aabb3f test_aabb = new aabb3f(); - - DAABBC aabb = BoxGeom.getAABB(); - test_aabb.minX = (float) aabb.getMin0(); - test_aabb.maxX = (float) aabb.getMax0(); - test_aabb.minY = (float) aabb.getMin1(); - test_aabb.maxY = (float) aabb.getMax1(); - test_aabb.minZ = (float) aabb.getMin2(); - test_aabb.maxZ = (float) aabb.getMax2(); - - GimDynArrayInt collision_result = GimDynArrayInt.GIM_CREATE_BOXQUERY_LIST(); - - ptrimesh.getAabbSet().gim_aabbset_box_collision(test_aabb, collision_result); - - if(collision_result.size()==0) - { - collision_result.GIM_DYNARRAY_DESTROY(); - return 0; - } - //*****Set globals for box collision******// - - //collide triangles - - //GUINT32 * boxesresult = GIM_DYNARRAY_POINTER(GUINT32,collision_result); - int[] boxesresult = collision_result.GIM_DYNARRAY_POINTER(); - ptrimesh.gim_trimesh_locks_work_data(); - - int ctContacts0 = 0; - - DVector3[] dv = { new DVector3(), new DVector3(), new DVector3() }; - //vec3f[] dv = new vec3f[] { new vec3f(), new vec3f(), new vec3f() };//[3]; - for(int i=0;i maxcontacts) diff --git a/core/src/main/java/org/ode4j/ode/internal/CollisionLibccd.java b/core/src/main/java/org/ode4j/ode/internal/CollisionLibccd.java index 15844f09..84af0c36 100644 --- a/core/src/main/java/org/ode4j/ode/internal/CollisionLibccd.java +++ b/core/src/main/java/org/ode4j/ode/internal/CollisionLibccd.java @@ -42,6 +42,8 @@ import static org.ode4j.ode.internal.CommonEnums.*; import static org.ode4j.ode.internal.DxCollisionUtil.dQuatTransform; import static org.ode4j.ode.internal.libccd.CCD.*; +import static org.ode4j.ode.internal.libccd.CCDCustomQuat.ccdQuatRotVec2; +import static org.ode4j.ode.internal.libccd.CCDCustomVec3.*; import static org.ode4j.ode.internal.libccd.CCDMPR.*; import static org.ode4j.ode.internal.libccd.CCDQuat.*; import static org.ode4j.ode.internal.libccd.CCDVec3.*; @@ -197,9 +199,8 @@ static void ccdGeomToCyl(final DxCylinder g, ccd_cyl_t cyl) { ccdQuatRotVec(cyl.axis, cyl.rot); ccdVec3Copy(cyl.p1, cyl.axis); ccdVec3Copy(cyl.p2, cyl.axis); - boolean cylAxisNormalizationResult = ccdVec3Normalize(cyl.axis); - dUVERIFY(!cylAxisNormalizationResult, "Invalid cylinder has been passed"); - ccdVec3Normalize(cyl.axis); + int cylAxisNormalizationResult = ccdVec3SafeNormalize(cyl.axis); + dUVERIFY(cylAxisNormalizationResult == 0, "Invalid cylinder has been passed"); ccdVec3Scale(cyl.p2, -1.0); ccdVec3Add(cyl.p1, cyl.pos); ccdVec3Add(cyl.p2, cyl.pos); @@ -696,10 +697,9 @@ static boolean correctTriangleContactNormal(ccd_triangle_t t, DContactGeom conta // Triangle face normal ccd_vec3_t triNormal = new ccd_vec3_t(); ccdVec3Cross(triNormal, edges[dMTV_FIRST], edges[dMTV_SECOND]); - // if (ccdVec3Normalize(triNormal)) { - // anyFault = true; - // } - anyFault |= ccdVec3Normalize(triNormal); + if (ccdVec3SafeNormalize(triNormal) != 0) { + anyFault = true; + } // Check the edges to see if one of them is involved for (int testEdgeIndex = !anyFault ? dMTV__MIN : dMTV__MAX; testEdgeIndex != dMTV__MAX; ++testEdgeIndex) { @@ -707,7 +707,7 @@ static boolean correctTriangleContactNormal(ccd_triangle_t t, DContactGeom conta ccd_vec3_t edgeAxis = edges[testEdgeIndex]; // Edge axis - if (ccdVec3Normalize(edgeAxis)) { + if (ccdVec3SafeNormalize(edgeAxis) != 0) { // This should not happen normally as in the case on of edges is degenerated // the triangle normal calculation would have to fail above. If for some // reason the above calculation succeeds and this one would not, it is diff --git a/core/src/main/java/org/ode4j/ode/internal/Common.java b/core/src/main/java/org/ode4j/ode/internal/Common.java index 991b7028..f8caae3a 100644 --- a/core/src/main/java/org/ode4j/ode/internal/Common.java +++ b/core/src/main/java/org/ode4j/ode/internal/Common.java @@ -24,9 +24,6 @@ *************************************************************************/ package org.ode4j.ode.internal; -import java.io.PrintWriter; -import java.io.StringWriter; - import org.ode4j.math.DVector3; import org.ode4j.ode.DGeom.DNearCallback; import org.ode4j.ode.DWorld; @@ -72,7 +69,7 @@ public class Common extends OdeConstants { throw new RuntimeException("dDOUBLE == dSINGLE"); } if (dDOUBLE) { - //TODO use MIN_VALUE instead? IEEE 754 ... + // dEpsilon = FLT_EPSILON or DBL_EPSILON dEpsilon = Double.MIN_NORMAL; MAX_FLOAT = Double.MAX_VALUE; } else { @@ -99,8 +96,7 @@ public class Common extends OdeConstants { public static final int d_MEMORY_OUT_OF_MEMORY = 1; -//From config-defaults.h - //TODO ??? + //From config-defaults.h /** @deprecated TZ this can be removed? */ @Deprecated public static final boolean dATOMICS_ENABLED = false; @@ -108,8 +104,6 @@ public class Common extends OdeConstants { public static final boolean dTRIMESH_OPCODE_USE_OLD_TRIMESH_TRIMESH_COLLIDER = false; - //TODO - //http://www.codeguru.com/forum/printthread.php?t=323835 //TODO use MACRO //#define EPSILON 0.0001 // Define your own tolerance //#define FLOAT_EQ(x,v) (((v - EPSILON) < x) && (x <( v + EPSILON))) @@ -161,11 +155,7 @@ public static void dUASSERT(int a, String msg) { } public static void dDEBUGMSG(String msg) { - StringWriter sw = new StringWriter(); - new PrintWriter(sw); - new RuntimeException(msg).printStackTrace(new PrintWriter(sw)); - String msg2 = sw.toString(); - dMessage (d_ERR_UASSERT, msg2); + dMessage (d_ERR_UASSERT, msg); } // #ifdef __GNUC__ @@ -185,7 +175,6 @@ public static void dDEBUGMSG(String msg) { // #endif public static void dSASSERT(boolean a) { - assert(a); if (!a) { dDebug (d_ERR_SASSERT, "Static assert failed"); } @@ -302,7 +291,7 @@ public static void dAVERIFY(Object a) { * @return Padded offset */ public static int dPAD(int a) { - return (((a) > 1) ? (((a) + 3) & (int)(~3)) : (a)); + return (((a) > 1) ? (((a) + 3) & (~3)) : (a)); } //#define dPAD(a) (((a) > 1) ? (((a) + 3) & (int)(~3)) : (a)) @@ -498,21 +487,21 @@ static class dDynamicsAxis { //#elif defined(dDOUBLE) //#define REAL(x) (x) - //#define dRecip(x) (1.0/(x)) //TODO replace + //#define dRecip(x) (1.0/(x)) public static double dRecip(double x) { return 1.0/x; } - //#define dSqrt(x) sqrt(x) //TODO replace + //#define dSqrt(x) sqrt(x) public static double dSqrt(double x) { return Math.sqrt(x); } //#define dRecipSqrt(x) (1.0/sqrt(x)) public static double dRecipSqrt(double x) { return 1.0/Math.sqrt(x); } - //#define dSin(x) sin(x)//TODO replace + //#define dSin(x) sin(x) public static double dSin(double x) { return Math.sin(x); } - //#define dCos(x) cos(x) //TODO replace + //#define dCos(x) cos(x) public static double dCos(double x) { return Math.cos(x); } - //#define dFabs(x) fabs(x) //TODO replace + //#define dFabs(x) fabs(x) public static double dFabs(double x) { return Math.abs(x); } - //#define dAtan2(y,x) atan2((y),(x)) //TODO replace + //#define dAtan2(y,x) atan2((y),(x)) public static double dAtan2(double y, double x) { return Math.atan2(y, x); } @@ -525,12 +514,12 @@ public static double dAcos(double x) { return Math.acos(x); } - //#define dFMod(a,b) (fmod((a),(b))) //TODO replace + //#define dFMod(a,b) (fmod((a),(b))) public static double dFMod(double x) { throw new UnsupportedOperationException(); //return Math.fmod(x); } - //#define dFloor(x) floor(x) //TODO replace + //#define dFloor(x) floor(x) public static double dFloor(double x) { return Math.floor(x); } //#define dCeil(x) ceilf(x) /* ceil */ @@ -680,23 +669,21 @@ enum { public static int SIZE_MAX = Integer.MAX_VALUE; -//#ifndef offsetof -//#define offsetof(s, m) ((size_t)&(((s *)8)->m) - (size_t)8) -// #endif -//#ifndef membersize -//#define membersize(s, m) (sizeof(((s *)8)->m)) -// #endif -//#ifndef endoffsetof -//#define endoffsetof(s, m) ((size_t)((size_t)&(((s *)8)->m) - (size_t)8) + sizeof(((s *)8)->m)) -// #endif - -//#define dMACRO_MAX(a, b) ((a) > (b) ? (a) : (b)) -// #define dMACRO_MIN(a, b) ((a) < (b) ? (a) : (b)) - public static int dMACRO_MIN(int a, int b) { - return Math.min(a, b); - } -// -// #define dMAKE_PADDING_SIZE(DataType, ElementType) ((sizeof(DataType) + sizeof(ElementType) - 1) / sizeof(ElementType)) + //#ifndef offsetof + //#define offsetof(s, m) ((size_t)&(((s *)8)->m) - (size_t)8) + // #endif + //#ifndef membersize + //#define membersize(s, m) (sizeof(((s *)8)->m)) + // #endif + //#ifndef endoffsetof + //#define endoffsetof(s, m) ((size_t)((size_t)&(((s *)8)->m) - (size_t)8) + sizeof(((s *)8)->m)) + // #endif + + //#define dMACRO_MAX(a, b) ((a) > (b) ? (a) : (b)) + // #define dMACRO_MIN(a, b) ((a) < (b) ? (a) : (b)) + // public static int dMACRO_MIN(int a, int b) { + // return Math.min(a, b); + // } // template @@ -766,10 +753,10 @@ public static void dxSwap(boolean[] one, int oneP, boolean[] another, int anothe // return value < lo ? (value_type)lo : value > hi ? (value_type)hi : value; // } static double dxClamp(double value, double lo, double hi) { - return value < lo ? lo : value > hi ? hi : value; + return value < lo ? lo : Math.min(value, hi); } static int dxClamp(int value, int lo, int hi) { - return value < lo ? lo : value > hi ? hi : value; + return value < lo ? lo : Math.min(value, hi); } diff --git a/core/src/main/java/org/ode4j/ode/internal/DLCP.java b/core/src/main/java/org/ode4j/ode/internal/DLCP.java index 8635bb81..0ce8fa53 100644 --- a/core/src/main/java/org/ode4j/ode/internal/DLCP.java +++ b/core/src/main/java/org/ode4j/ode/internal/DLCP.java @@ -29,7 +29,10 @@ import static org.ode4j.ode.internal.Common.*; import static org.ode4j.ode.internal.ErrorHandler.dDebug; import static org.ode4j.ode.internal.ErrorHandler.dMessage; -import static org.ode4j.ode.internal.FastLDLT.dxtFactorLDLT; +import static org.ode4j.ode.internal.FastLDLTFactor.factorMatrixAsLDLT; +import static org.ode4j.ode.internal.FastLDLTSolve.solveEquationSystemWithLDLT; +import static org.ode4j.ode.internal.FastLSolve.solveL1Straight; +import static org.ode4j.ode.internal.FastLTSolve.solveL1Transposed; import static org.ode4j.ode.internal.Matrix.*; import static org.ode4j.ode.internal.Misc.dClearUpperTriangle; import static org.ode4j.ode.internal.Misc.dMakeRandomMatrix; @@ -205,7 +208,7 @@ public class DLCP { //#define DEBUG_LCP //#define dLCP_FAST // use fast dLCP object - //private static final boolean dLCP_FAST = true; + private static final boolean dLCP_FAST = true; //#define NUB_OPTIMIZATIONS protected static final boolean NUB_OPTIMIZATIONS = true; @@ -262,9 +265,196 @@ static void transfer_b_to_x(double[] pairsbxA, int pairsbxP, int n, boolean zero } } + /** + * swap row/column i1 with i2 in the n*n matrix A. the leading dimension of + * A is nskip. this only references and swaps the lower triangle. + * if `do_fast_row_swaps' is nonzero and row pointers are being used, then + * rows will be swapped by exchanging row pointers. otherwise the data will + * be copied. + */ + private static void swapRowsAndCols (double[] A, int n, int i1, int i2, int nskip, + boolean do_fast_row_swaps) + { + //dAASSERT (A); + dAASSERT (n > 0 && i1 >= 0 && i2 >= 0 && i1 < n && i2 < n && + nskip >= n && i1 < i2); + + if (ROWPTRS) {//# ifdef ROWPTRS +// dReal *A_i1 = A[i1]; +// dReal *A_i2 = A[i2]; +// for (int i=i1+1; i0 && i1 < n && i2 < n && nskip >= n && i1 <= i2); + + if (i1 != i2) { + swapRowsAndCols (A, n, i1, i2, nskip, do_fast_row_swaps); + +// dxSwap((pairsbx + (size_t)i1 * PBX__MAX)[PBX_B], (pairsbx + (size_t)i2 * PBX__MAX)[PBX_B]); +// dxSwap((pairsbx + (size_t)i1 * PBX__MAX)[PBX_X], (pairsbx + (size_t)i2 * PBX__MAX)[PBX_X]); +// dSASSERT(PBX__MAX == 2); + dxSwap(pairsbxA, pairsbxP + i1 * PBX__MAX + PBX_B, pairsbxA, pairsbxP + i2 * PBX__MAX + PBX_B); + dxSwap(pairsbxA, pairsbxP + i1 * PBX__MAX + PBX_X, pairsbxA, pairsbxP + i2 * PBX__MAX + PBX_X); + dSASSERT(PBX__MAX == 2); + + dxSwap(w, i1, w, i2); + +// dxSwap((pairslh + (size_t)i1 * PLH__MAX)[PLH_LO], (pairslh + (size_t)i2 * PLH__MAX)[PLH_LO]); +// dxSwap((pairslh + (size_t)i1 * PLH__MAX)[PLH_HI], (pairslh + (size_t)i2 * PLH__MAX)[PLH_HI]); +// dSASSERT(PLH__MAX == 2); + dxSwap(pairslhA, pairslhP + i1 * PLH__MAX + PLH_LO, pairslhA, pairslhP + i2 * PLH__MAX + PLH_LO); + dxSwap(pairslhA, pairslhP + i1 * PLH__MAX + PLH_HI, pairslhA, pairslhP + i2 * PLH__MAX + PLH_HI); + dSASSERT(PLH__MAX == 2); + + dxSwap(p, i1, p, i2); + dxSwap(state, i1, state, i2); + + if (findex != null) { + dxSwap(findex, i1, findex, i2); + } + } + } + + + // for debugging - check that L,d is the factorization of A[C,C]. + // A[C,C] has size nC*nC and leading dimension nskip. + // L has size nC*nC and leading dimension nskip. + // d has size nC. + + //#ifdef DEBUG_LCP -> TZ see first 'return' + + //static void checkFactorization (ATYPE A, dReal *_L, dReal *_d, + // int nC, int *C, int nskip) + protected void checkFactorization (double[] A, double[]_L, double[]_d, + int nC, int[] C, int nskip) + { + if (!DEBUG_LCP) return; + int i,j; + if (nC==0) return; + + // get A1=A, copy the lower triangle to the upper triangle, get A2=A[C,C] + DMatrixN A1 = new DMatrixN(nC,nC); + for (i=0; i 1e-8) + dDebug (0,"L*D*L' check, maximum difference = %.6e\n",diff); + } + + //#endif + + + // for debugging + + //#ifdef DEBUG_LCP -> TZ see first 'return' + + //static void checkPermutations (int i, int n, int nC, int nN, int *p, int *C) + protected static void checkPermutations (int i, int n, int nC, int nN, + int[] p, int[] C) + { + if (!DEBUG_LCP) return; + int j,k; + dIASSERT(/*nC>=0 && nN>=0 &&*/ (nC + nN) == i && i < n); + for (k=0; k= 0 && p[k] < i); + for (k=i; k(m_pairsbx + PBX_X, m_n); dxtSetZero(m_pairsbxA, 0 + PBX_X, m_n, PBX__MAX); @@ -493,8 +684,8 @@ static int estimate_transfer_i_from_C_to_N_mem_req(int nC, int nskip) { } } transfer_b_to_x(m_pairsbxA, 0, nub, false); - dxtFactorLDLT(m_L, m_d, nub, m_nskip, 1); - dxtSolveLDLT(m_L, m_d, 0, m_pairsbxA, 0 + PBX_X, nub, m_nskip, 1, PBX__MAX); + factorMatrixAsLDLT(m_L, m_d, nub, m_nskip, 1); + solveEquationSystemWithLDLT(m_L, m_d, 0, m_pairsbxA, 0 + PBX_X, nub, m_nskip, 1, PBX__MAX); dSetZero (m_w,nub); { int[] C = m_C; @@ -585,7 +776,7 @@ void transfer_i_from_N_to_C (int i) for (int j=0; j= lo[i] && xi <= hi[i] && wi == 0) { return 1; } - - /** - * swap row/column i1 with i2 in the n*n matrix A. the leading dimension of - * A is nskip. this only references and swaps the lower triangle. - * if `do_fast_row_swaps' is nonzero and row pointers are being used, then - * rows will be swapped by exchanging row pointers. otherwise the data will - * be copied. - */ - private static void swapRowsAndCols (double[] A, int n, int i1, int i2, int nskip, - boolean do_fast_row_swaps) - { - //dAASSERT (A); - dAASSERT (n > 0 && i1 >= 0 && i2 >= 0 && i1 < n && i2 < n && - nskip >= n && i1 < i2); - - if (ROWPTRS) {//# ifdef ROWPTRS -// dReal *A_i1 = A[i1]; -// dReal *A_i2 = A[i2]; -// for (int i=i1+1; i0 && i1 < n && i2 < n && nskip >= n && i1 <= i2); - - if (i1 != i2) { - swapRowsAndCols (A, n, i1, i2, nskip, do_fast_row_swaps); - -// dxSwap((pairsbx + (size_t)i1 * PBX__MAX)[PBX_B], (pairsbx + (size_t)i2 * PBX__MAX)[PBX_B]); -// dxSwap((pairsbx + (size_t)i1 * PBX__MAX)[PBX_X], (pairsbx + (size_t)i2 * PBX__MAX)[PBX_X]); -// dSASSERT(PBX__MAX == 2); - dxSwap(pairsbxA, pairsbxP + i1 * PBX__MAX + PBX_B, pairsbxA, pairsbxP + i2 * PBX__MAX + PBX_B); - dxSwap(pairsbxA, pairsbxP + i1 * PBX__MAX + PBX_X, pairsbxA, pairsbxP + i2 * PBX__MAX + PBX_X); - dSASSERT(PBX__MAX == 2); - - dxSwap(w, i1, w, i2); - -// dxSwap((pairslh + (size_t)i1 * PLH__MAX)[PLH_LO], (pairslh + (size_t)i2 * PLH__MAX)[PLH_LO]); -// dxSwap((pairslh + (size_t)i1 * PLH__MAX)[PLH_HI], (pairslh + (size_t)i2 * PLH__MAX)[PLH_HI]); -// dSASSERT(PLH__MAX == 2); - dxSwap(pairslhA, pairslhP + i1 * PLH__MAX + PLH_LO, pairslhA, pairslhP + i2 * PLH__MAX + PLH_LO); - dxSwap(pairslhA, pairslhP + i1 * PLH__MAX + PLH_HI, pairslhA, pairslhP + i2 * PLH__MAX + PLH_HI); - dSASSERT(PLH__MAX == 2); - - dxSwap(p, i1, p, i2); - dxSwap(state, i1, state, i2); - - if (findex != null) { - dxSwap(findex, i1, findex, i2); - } - } - } - - - // for debugging - check that L,d is the factorization of A[C,C]. - // A[C,C] has size nC*nC and leading dimension nskip. - // L has size nC*nC and leading dimension nskip. - // d has size nC. - - //#ifdef DEBUG_LCP -> TZ see first 'return' - - //static void checkFactorization (ATYPE A, dReal *_L, dReal *_d, - // int nC, int *C, int nskip) - protected void checkFactorization (double[] A, double[]_L, double[]_d, - int nC, int[] C, int nskip) - { - if (!DEBUG_LCP) return; - int i,j; - if (nC==0) return; - - // get A1=A, copy the lower triangle to the upper triangle, get A2=A[C,C] - DMatrixN A1 = new DMatrixN(nC,nC); - for (i=0; i 1e-8) - dDebug (0,"L*D*L' check, maximum difference = %.6e\n",diff); - } - - //#endif - - - // for debugging - - //#ifdef DEBUG_LCP -> TZ see first 'return' - - //static void checkPermutations (int i, int n, int nC, int nN, int *p, int *C) - protected static void checkPermutations (int i, int n, int nC, int nN, - int[] p, int[] C) - { - if (!DEBUG_LCP) return; - int j,k; - dIASSERT(/*nC>=0 && nN>=0 &&*/ (nC + nN) == i && i < n); - for (k=0; k= 0 && p[k] < i); - for (k=i; k0) - { - out = true; - break; - }; + out = true; + } + + if (!out) { + out = false; + for (int k = 0; k < cvx1.planecount; ++k) { + // d = p[0]*cvx1.planes[(k*4)+0]+ + // p[1]*cvx1.planes[(k*4)+1]+ + // p[2]*cvx1.planes[(k*4)+2]- + // cvx1.planes[(k*4)+3]; + d = p.dot(cvx1.planesV[k]) - cvx1.planesD[k]; + if (d > 0) { + out = true; + break; + } + ; + } } if(!out) { @@ -1632,17 +1630,13 @@ else if(!CheckSATConvexEdges(cvx1,cvx2,ccso)) d = p.dot(rplaneV) - rplaneD; if(d>0) { - //dVector3Copy(p,SAFECONTACT(flags, contact, contacts, skip).pos); - //dVector3Copy(rplane,SAFECONTACT(flags, contact, contacts, skip).normal); - //SAFECONTACT(flags, contact, contacts, skip).g1=cvx1;//&cvx1; - //SAFECONTACT(flags, contact, contacts, skip).g2=cvx2;//&cvx2; - //SAFECONTACT(flags, contact, contacts, skip).depth=d; - DContactGeom contact = contactBuf.getSafe(flags, contacts); - contact.pos.set(p); - contact.normal.set(rplaneV); - contact.g1 = cvx1; - contact.g2 = cvx2; - contact.depth = d; + // SAFECONTACT(flags, contact, contacts, skip); + DContactGeom target = contactBuf.getSafe(flags, contacts); + dVector3Copy(p, target.pos); + dVector3Copy(rplaneV, target.normal); + target.g1 = cvx1; + target.g2 = cvx2; + target.depth = d; ++contacts; if (contacts==maxc) return contacts; } @@ -1661,17 +1655,13 @@ else if(!CheckSATConvexEdges(cvx1,cvx2,ccso)) d = i1.dot(rplaneV) - rplaneD; if(d>0) { - DContactGeom contact = contactBuf.getSafe(flags, contacts); - //dVector3Copy(i1,SAFECONTACT(flags, contact, contacts, skip).pos); - //dVector3Copy(rplane,SAFECONTACT(flags, contact, contacts, skip).normal); - //SAFECONTACT(flags, contact, contacts, skip).g1=cvx1;//&cvx1; - //SAFECONTACT(flags, contact, contacts, skip).g2=cvx2;//&cvx2; - //SAFECONTACT(flags, contact, contacts, skip).depth=d; - contact.pos.set(i1); - contact.normal.set(rplaneV); - contact.g1 = cvx1; - contact.g2 = cvx2; - contact.depth = d; + // dContactGeom *target = SAFECONTACT(flags, contact, contacts, skip); + DContactGeom target = contactBuf.getSafe(flags, contacts); + dVector3Copy(i1, target.pos); + dVector3Copy(rplaneV, target.normal); + target.g1 = cvx1; + target.g2 = cvx2; + target.depth = d; ++contacts; if (contacts==maxc) return contacts; } @@ -1739,10 +1729,9 @@ else if(!CheckSATConvexEdges(cvx1,cvx2,ccso)) outside = false; for(int j=0;j0) { - //dVector3Copy(i1,SAFECONTACT(flags, contact, contacts, skip).pos); - //dVector3Copy(rplane,SAFECONTACT(flags, contact, contacts, skip).normal); - //SAFECONTACT(flags, contact, contacts, skip).g1=cvx1;//&cvx1; - //SAFECONTACT(flags, contact, contacts, skip).g2=cvx2;//&cvx2; - //SAFECONTACT(flags, contact, contacts, skip).depth=d; - DContactGeom contact = contactBuf.getSafe(flags, contacts); - contact.pos.set(i1); - contact.normal.set(rplaneV); - contact.g1 = cvx1; - contact.g2 = cvx2; - contact.depth = d; + // dContactGeom *target = SAFECONTACT(flags, contact, contacts, skip); + DContactGeom target = contactBuf.getSafe(flags, contacts); + dVector3Copy(i1, target.pos); + dVector3Copy(rplaneV, target.normal); + target.g1 = cvx1; + target.g2 = cvx2; + target.depth = d; ++contacts; if (contacts==maxc) return contacts; } @@ -1774,51 +1759,38 @@ else if(!CheckSATConvexEdges(cvx1,cvx2,ccso)) } } } - else if(ccso.depth_type==2) // edge-edge + else if (ccso.depth_type == 2) // edge-edge { -// // Some parts borrowed from dBoxBox -// DVector3 ua = new DVector3(),ub = new DVector3(),pa = new DVector3(),pb = new DVector3(); -// RefDouble alpha=new RefDouble(),beta=new RefDouble(); -// // Get direction of first edge -// //for (i=0; i<3; i++) ua[i] = ccso.e1b[i]-ccso.e1a[i]; -// ua.eqDiff(ccso.e1b, ccso.e1a); -// dNormalize3(ua); // normalization shouldn't be necesary but dLineClosestApproach requires it -// // Get direction of second edge -// //for (i=0; i<3; i++) ub[i] = ccso.e2b[i]-ccso.e2a[i]; -// ub.eqDiff(ccso.e2b, ccso.e2a); -// dNormalize3(ub); // same as with ua normalization -// // Get closest points between edges (one at each) -// DxCollisionUtil.dLineClosestApproach (ccso.e1a,ua,ccso.e2a,ub,alpha,beta); -// //for (i=0; i<3; i++) pa[i] = ccso.e1a[i]+(ua[i]*alpha.get()); -// pa.eqSum(ccso.e1a, ua, alpha.get()); -// //for (i=0; i<3; i++) pb[i] = ccso.e2a[i]+(ub[i]*beta.get()); -// pb.eqSum(ccso.e2a, ub, beta.get()); -// // Set the contact point as halfway between the 2 closest points -//// for (i=0; i<3; i++) SAFECONTACT(flags, contact, contacts, skip).pos[i] = REAL(0.5)*(pa[i]+pb[i]); -//// SAFECONTACT(flags, contact, contacts, skip).g1=cvx1;//&cvx1; -//// SAFECONTACT(flags, contact, contacts, skip).g2=cvx2;//&cvx2; -//// dVector3Copy(ccso.plane,SAFECONTACT(flags, contact, contacts, skip).normal); -//// SAFECONTACT(flags, contact, contacts, skip).depth=ccso.min_depth; -// DContactGeom contact = contactBuf.getSafe(flags, contacts); -// contact.pos.eqSum(pa, 0.5, pb, 0.5); -// contact.g1 = cvx1; -// contact.g2 = cvx2; -// //TODO TZ optimize with dVector3! -// contact.normal.set(ccso.plane.get0(), ccso.plane.get1(), ccso.plane.get2()); -// contact.depth = ccso.min_depth; -// ++contacts; - DVector3 c1 = new DVector3(), c2 = new DVector3(); - DContactGeom contact = contactBuf.getSafe(flags, contacts); - //SAFECONTACT(flags, contact, contacts, skip)->depth = - contact.depth = dSqrt(ClosestPointBetweenSegments(ccso.e1a,ccso.e1b,ccso.e2a,ccso.e2b,c1,c2)); - contact.g1=cvx1; - contact.g2=cvx2; - //dVector3Copy(c1,SAFECONTACT(flags, contact, contacts, skip)->pos); - contact.pos.set( c1 ); - contact.normal.eqDiff(c2, c1); - contact.normal.normalize(); - //dNormalize3(SAFECONTACT(flags, contact, contacts, skip)->normal); + ClosestPointBetweenSegments(ccso.e1a, ccso.e1b, ccso.e2a, ccso.e2b, c1, c2); + + // dContactGeom *target = SAFECONTACT(flags, contact, contacts, skip); + DContactGeom target = contactBuf.getSafe(flags, contacts); + dSubtractVectors3(target.normal, c2, c1); + double depth_square = dCalcVectorLengthSquare3(target.normal); + + if (dxSafeNormalize3(target.normal)) + { + target.depth = dSqrt(depth_square); + } + else + { + // If edges coincide return direction from one center to the other as the contact normal + dVector3Copy(ccso.dist, target.normal); + + if (!dxSafeNormalize3(target.normal)) + { + // If the both centers coincide as well return an arbitrary vector. The depth is going to be zero anyway. + // dAssignVector3(target.normal, 1, 0, 0); + target.normal.set(1, 0, 0); + } + + target.depth = 0; // Since the edges coincide, return a contact of zero depth + } + + target.g1 = cvx1; + target.g2 = cvx2; + dVector3Copy(c1, target.pos); contacts++; } return contacts; @@ -1936,17 +1908,31 @@ int dCollideRayConvex( DxRay ray, DxConvex convex, contact.side2 = -1; // TODO: set plane index? double alpha, beta, nsign; + boolean flag = false; // // Compute some useful info // - - DVector3 ray_pos = new DVector3(); - DVector3 ray_dir = new DVector3(); - dMultiply1_331(ray_pos, convex.final_posr().R(), new DVector3(ray.final_posr().pos()).sub(convex.final_posr().pos())); - dMultiply1_331(ray_dir, convex.final_posr().R(), ray.final_posr().R().columnAsNewVector(2)); - boolean flag = false; + // dVector3 ray_pos = { + // ray->final_posr->pos[0] - convex->final_posr->pos[0], + // ray->final_posr->pos[1] - convex->final_posr->pos[1], + // ray->final_posr->pos[2] - convex->final_posr->pos[2] + // }; + DVector3 ray_pos = new DVector3(ray.final_posr().pos()).sub(convex.final_posr().pos()); + // + // dVector3 ray_dir = { + // ray->final_posr->R[0 * 4 + 2], + // ray->final_posr->R[1 * 4 + 2], + // ray->final_posr->R[2 * 4 + 2] + // }; + DVector3 ray_dir = ray.final_posr().R().columnAsNewVector(2); + + //dMultiply1_331(ray_pos, convex->final_posr->R, ray_pos); + //dMultiply1_331(ray_dir, convex->final_posr->R, ray_dir); + dMultiply1_331(ray_pos, convex.final_posr().R(), ray_pos); + dMultiply1_331(ray_dir, convex.final_posr().R(), ray_dir); + for ( int i = 0; i < convex.planecount; ++i ) { // Alias this plane. diff --git a/core/src/main/java/org/ode4j/ode/internal/DxHashSpace.java b/core/src/main/java/org/ode4j/ode/internal/DxHashSpace.java index e9cf874e..5c524cde 100644 --- a/core/src/main/java/org/ode4j/ode/internal/DxHashSpace.java +++ b/core/src/main/java/org/ode4j/ode/internal/DxHashSpace.java @@ -226,12 +226,15 @@ public void collide (Object data, DNearCallback callback) aabb.level = level; if (level > maxlevel) maxlevel = level; // cellsize = 2^level - double cellSizeRecip = dRecip(ldexp(1.0, level)); - // discretize AABB position to cell size + double cellSizeRecip = dRecip(ldexp(1.0, level)); // No computational errors here! // discretize AABB position to cell size for (i=0; i < 3; i++) { - aabb.dbounds[2*i] = (int)Math.floor (geom._aabb.getMin(i) * cellSizeRecip); - aabb.dbounds[2*i+1] = (int)Math.floor (geom._aabb.getMax(i) * cellSizeRecip); + double aabbBoundMin = Math.floor (geom._aabb.getMin(i) * cellSizeRecip); // No computational errors so far! + double aabbBoundMax = Math.floor (geom._aabb.getMax(i) * cellSizeRecip); // No computational errors so far! + dICHECK(aabbBoundMin >= Integer.MIN_VALUE && aabbBoundMin = Integer.MIN_VALUE && aabbBoundMax = RANDOM_CONSTRAINTS_REORDERING_FREQUENCY + ? (iteration % RANDOM_CONSTRAINTS_REORDERING_FREQUENCY < RRS_MAX) + : iteration == 0) { result = true; } } else { @@ -911,22 +936,19 @@ void dxQuickStepIsland_Stage1(dxQuickStepperStage1CallContext stage1CallContext) { int mcurrO = 0;//mindex; int moffs = 0, mfboffs = 0; - mindex[mcurrO+0] = moffs; - mindex[mcurrO+1] = mfboffs; - mcurrO += 2; + setMIndex(mindex, mcurrO, moffs, mfboffs); + ++mcurrO; + //for (DJointWithInfo1 jicurr: jointinfos) { for (int i = 0; i < nj; i++) { DJointWithInfo1 jicurr = jointinfos[i]; - //TODO fix issue #18 - //for (DJointWithInfo1 jicurr: jointinfos) { //const dJointWithInfo1 *const jiend = jointinfos + nj; //for (const dJointWithInfo1 *jicurr = jointinfos; jicurr != jiend; ++jicurr) { DxJoint joint = jicurr.joint; moffs += jicurr.info.m; if (joint.feedback != null) { mfboffs += jicurr.info.m; } - mindex[mcurrO+0] = moffs; - mindex[mcurrO+1] = mfboffs; - mcurrO += 2; + setMIndex(mindex, mcurrO, moffs, mfboffs); + ++mcurrO; } } @@ -1095,19 +1117,21 @@ void dxQuickStepIsland_Stage2a(dxQuickStepperStage2CallContext stage2CallContext int ji; while ((ji = Atomics.ThrsafeIncrementIntUpToLimit(stage2CallContext.m_ji_J, nj)) != nj) { - final int ofsi = mindex[ji * 2 + 0]; - final int infom = mindex[ji * 2 + 2] - ofsi; + final int ofsi = getMIndex(mindex, ji); + final int infom = getMIndex(mindex, ji + 1) - ofsi; int jRow = ofsi * JME__MAX; - int jEnd = jRow + infom * JME__MAX; - for (int jCurr = jRow; jCurr != jEnd; jCurr += JME__MAX) { - dSetZero(J, jCurr + JME__J1_MIN, JME__J1_COUNT); - J[jCurr + JME_RHS] = 0.0; - J[jCurr + JME_CFM] = worldCFM; - dSetZero(J, jCurr + JME__J2_MIN, JME__J2_COUNT); - J[jCurr + JME_LO] = -dInfinity; - J[jCurr + JME_HI] = dInfinity; - } - dSASSERT(JME__J1_COUNT + 2 + JME__J2_COUNT + 2 == JME__MAX); + { + int jEnd = jRow + infom * JME__MAX; + for (int jCurr = jRow; jCurr != jEnd; jCurr += JME__MAX) { + dSetZero(J, jCurr + JME__J1_MIN, JME__J1_COUNT); + J[jCurr + JME_RHS] = 0.0; + J[jCurr + JME_CFM] = worldCFM; + dSetZero(J, jCurr + JME__J2_MIN, JME__J2_COUNT); + J[jCurr + JME_LO] = -dInfinity; + J[jCurr + JME_HI] = dInfinity; + } + dSASSERT(JME__J1_COUNT + 2 + JME__J2_COUNT + 2 == JME__MAX); + } int findexRow = ofsi; dSetValue(findex, findexRow, infom, -1); @@ -1116,41 +1140,72 @@ void dxQuickStepIsland_Stage2a(dxQuickStepperStage2CallContext stage2CallContext JME__MAX, J, jRow + JME__RHS_CFM_MIN, J, jRow + JME__LO_HI_MIN, findex, findexRow); // findex iteration is compact and is not going to pollute caches - do it first - // adjust returned findex values for global index numbering - // TZ: This looks alright, but is different from the original where we use findexRow - for (int j = infom; j != 0; ) { - --j; - int fival = findex[j+ofsi]; - if (fival != -1) { - findex[j+ofsi] = fival + ofsi; - validFIndices++; - } - } - - for (int jCurr = jRow; jCurr != jEnd; jCurr += JME__MAX) { - J[jCurr + JME_RHS] *= stepsizeRecip; - J[jCurr + JME_CFM] *= stepsizeRecip; - } - + { + // adjust returned findex values for global index numbering + int findicesEndOfs = findexRow + infom; + for (int findexCurrOfs = findexRow; findexCurrOfs != findicesEndOfs; ++findexCurrOfs) { + int fival = findex[findexCurrOfs]; // *findexCurr; + if (fival != -1) { + findex[findexCurrOfs] = fival + ofsi; // *findexCurr = fival + ofsi; + ++validFIndices; + } + } + // TZ: This looks alright, but is different from the original where we use findexRow +// for (int j = infom; j != 0; ) { +// --j; +// int fival = findex[j + ofsi]; +// if (fival != -1) { +// findex[j + ofsi] = fival + ofsi; +// ++validFIndices; +// } +// } + } + { + // dReal *const JEnd = JRow + infom * JME__MAX; + int jEnd = jRow + infom * JME__MAX; + for (int jCurr = jRow; jCurr != jEnd; jCurr += JME__MAX) { + J[jCurr + JME_RHS] *= stepsizeRecip; + J[jCurr + JME_CFM] *= stepsizeRecip; + } + } // we need a copy of Jacobian for joint feedbacks // because it gets destroyed by SOR solver // instead of saving all Jacobian, we can save just rows // for joints, that requested feedback (which is normally much less) - int mfbcurr = mindex[ji * 2 + 1], mfbnext = mindex[ji * 2 + 3]; - int mfbCount = mfbnext - mfbcurr; - if (mfbCount != 0) { - int jEndMfb = jRow + mfbCount * JME__MAX; - for (int jCurr = jRow; jCurr < jEndMfb ;jCurr += JME__MAX) { - System.arraycopy(J, jCurr + JME__J1_MIN, Jcopy, jCopy_ofs + JCE__J1_MIN, JME__J1_COUNT); - System.arraycopy(J, jCurr + JME__J2_MIN, Jcopy, jCopy_ofs + JCE__J2_MIN, JME__J2_COUNT); - jCopy_ofs += JCE__MAX; - } + int mfbIndex = getFbIndex(mindex, ji); //mindex[ji].fbIndex; + if (mfbIndex != getMIndex(mindex, ji + 1)) { //mindex[ji + 1].fbIndex) { + // dReal *const JEnd = JRow + infom * JME__MAX; + // dReal *JCopyRow = JCopy + mfbIndex * JCE__MAX; // Random access by mfbIndex here! Do not optimize! + final int jEnd = jRow + infom * JME__MAX; + int JCopyRow = jCopy_ofs + mfbIndex * JCE__MAX; // Random access by mfbIndex here! Do not optimize! + //for (const dReal *JCurr = JRow; ; ) { + for (int jCurr = jRow; jCurr < jEnd ;) { + //for (unsigned i = 0; i != JME__J1_COUNT; ++i) { JCopyRow[i + JCE__J1_MIN] = JCurr[i + JME__J1_MIN]; } + //for (unsigned j = 0; j != JME__J2_COUNT; ++j) { JCopyRow[j + JCE__J2_MIN] = JCurr[j + JME__J2_MIN]; } + System.arraycopy(J, jCurr + JME__J1_MIN, Jcopy, JCopyRow + JCE__J1_MIN, JME__J1_COUNT); + System.arraycopy(J, jCurr + JME__J2_MIN, Jcopy, JCopyRow + JCE__J2_MIN, JME__J2_COUNT); + JCopyRow += JCE__MAX; + + // TZ move this outside of loop + // dSASSERT((unsigned)JCE__J1_COUNT == JME__J1_COUNT); + // dSASSERT((unsigned)JCE__J2_COUNT == JME__J2_COUNT); + // dSASSERT(JCE__J1_COUNT + JCE__J2_COUNT == JCE__MAX); + + if ((jCurr += JME__MAX) == jEnd) { + break; + } + } } } Atomics.ThrsafeAdd(localContext.m_valid_findices, validFIndices); } - { + // TZ moved this here from inside loop above. + dSASSERT(JCE__J1_COUNT == JME__J1_COUNT); + dSASSERT(JCE__J2_COUNT == JME__J2_COUNT); + dSASSERT(JCE__J1_COUNT + JCE__J2_COUNT == JCE__MAX); + + { int[] jb = localContext.m_jb; // create an array of body numbers for each joint row @@ -1161,8 +1216,8 @@ void dxQuickStepIsland_Stage2a(dxQuickStepperStage2CallContext stage2CallContext int b1 = (joint.node[0].body!=null) ? (joint.node[0].body.tag) : -1; int b2 = (joint.node[1].body!=null) ? (joint.node[1].body.tag) : -1; - int jb_end = 2 * mindex[ji * 2 + 2]; - int jb_ptr = 2 * mindex[ji * 2 + 0]; + int jb_end = 2 * getMIndex(mindex, ji + 1); + int jb_ptr = 2 * getMIndex(mindex, ji); for (; jb_ptr != jb_end; jb_ptr += 2) { jb[jb_ptr] = b1; jb[jb_ptr+1] = b2; @@ -1407,8 +1462,8 @@ void dxQuickStepIsland_Stage4a(dxQuickStepperStage4CallContext stage4CallContext int ji_step; while ((ji_step = Atomics.ThrsafeIncrementIntUpToLimit(stage4CallContext.m_ji_4a, nj_steps)) != nj_steps) { int ji = ji_step * step_size; - int lambdacurr = mindex[ji * 2]; - int lambdsnext = mindex[2 * (ji + Math.min(step_size, nj - ji))]; + int lambdacurr = getMIndex(mindex, ji); + int lambdsnext = getMIndex(mindex, ji + Math.min(step_size, nj - ji)); dSetZero(lambda, lambdacurr, lambdsnext - lambdacurr); } } @@ -1602,11 +1657,12 @@ void dxQuickStepIsland_Stage4LCP_ReorderPrep(dxQuickStepperStage4CallContext sta { dxQuickStepperLocalContext localContext = stage4CallContext.m_localContext; int m = localContext.m_m; + int valid_findices = localContext.m_valid_findices.get(); IndexError []order = stage4CallContext.m_order; // make sure constraints with findex < 0 come first. int orderhead = 0; - int ordertail = m; + int ordertail = m - valid_findices; int[] findex = localContext.m_findex; // Fill the array from both ends @@ -1615,13 +1671,8 @@ void dxQuickStepIsland_Stage4LCP_ReorderPrep(dxQuickStepperStage4CallContext sta order[orderhead].index = i; // Place them at the front ++orderhead; } else { - // TODO CHECK-TZ Has this been fixed with the recent bug fixes? - // WARNING!!! - // The dependent constraints are put in reverse order here (backwards to front). - // They MUST be ordered this way. - // Putting them in normal order makes simulation less stable for some mysterious reason. - --ordertail; - order[ordertail].index = i; // Place them at the end + order[ordertail].index = i; // Place them at the end + ++ordertail; } } } @@ -1664,7 +1715,10 @@ int dxQuickStepIsland_Stage4LCP_IterationStart(final dxQuickStepperStage4CallCon reorderRequired = true; } - // Increment iterations counter in advance as anyway it needs to be incremented + // TODO CHECK-TZ why don't we need this? + // unsigned syncCallDependencies = reorderRequired ? 1 : stage4LCP_Iteration_allowedThreads; + + // Increment iterations counter in advance as anyway it needs to be incremented // before independent tasks (the reordering or the iteration) are posted // (otherwise next iteration may complete before the increment // and the same iteration index may be used again). @@ -1743,15 +1797,25 @@ void dxQuickStepIsland_Stage4LCP_ConstraintsReordering(dxQuickStepperStage4CallC dxQuickStepIsland_Stage4LCP_DependencyMapForNewOrderRebuilding(stage4CallContext); } } + else { + // TODO CHECK-TZ + if (true) + throw new UnsupportedOperationException("THis was added in 16.2, why did I miss it?"); + // NOTE: So far, this branch is only called in CONSTRAINTS_REORDERING_METHOD == REORDERING_METHOD__BY_ERROR case + if (Atomics.ThrsafeExchangeAdd(stage4CallContext.m_SOR_reorderThreadsRemaining, -1) == 1) { // If last thread has exited the reordering routine... + dIASSERT(iteration != 0); + dxQuickStepIsland_Stage4LCP_DependencyMapFromSavedLevelsReconstruction(stage4CallContext); + } + } } private static boolean dxQuickStepIsland_Stage4LCP_ConstraintsShuffling(dxQuickStepperStage4CallContext stage4CallContext, int iteration) { boolean result = false; + // #if CONSTRAINTS_REORDERING_METHOD == REORDERING_METHOD__BY_ERROR if (CONSTRAINTS_REORDERING_METHOD == ReorderingMethod.REORDERING_METHOD__BY_ERROR) { /* - #if CONSTRAINTS_REORDERING_METHOD == REORDERING_METHOD__BY_ERROR struct ConstraintsReorderingHelper { void operator ()(dxQuickStepperStage4CallContext *stage4CallContext, int int startIndex, int int endIndex) @@ -1784,6 +1848,13 @@ else if (last_lambda[i] != REAL(0.0)) { if (iteration > 1) { // Only reorder starting from iteration #2 // sort the constraints so that the ones converging slowest // get solved last. use the absolute (not relative) error. + // + // Full reorder needs to be done. + // Even though this contradicts the initial idea of moving dependent constraints + // to the order end the algorithm does not work the other way well. + // It looks like the iterative method needs a shake after it already found + // some initial approximations and those incurred errors help it to converge even better. + // if (ThrsafeExchange(&stage4CallContext->m_SOR_reorderHeadTaken, 1) == 0) { // Process the head const dxQuickStepperLocalContext *localContext = stage4CallContext->m_localContext; @@ -1822,37 +1893,55 @@ else if (iteration == 1) { else { result = true; // return true on 0th iteration to build dependency map for the initial order } - #elif CONSTRAINTS_REORDERING_METHOD == REORDERING_METHOD__RANDOMLY */ throw new UnsupportedOperationException(); + // #elif CONSTRAINTS_REORDERING_METHOD == REORDERING_METHOD__RANDOMLY } else if (CONSTRAINTS_REORDERING_METHOD == ReorderingMethod.REORDERING_METHOD__RANDOMLY) { if (iteration != 0) { - if (iteration % RANDOM_CONSTRAINTS_REORDERING_FREQUENCY == RRS_REORDERING) { + dIASSERT(!dIN_RANGE(iteration, 0, RANDOM_CONSTRAINTS_REORDERING_FREQUENCY)); + + dIASSERT(iteration % RANDOM_CONSTRAINTS_REORDERING_FREQUENCY == RRS_REORDERING); + { + dIASSERT(!dIN_RANGE(iteration, 0, RANDOM_CONSTRAINTS_REORDERING_FREQUENCY)); + + dIASSERT(iteration % RANDOM_CONSTRAINTS_REORDERING_FREQUENCY == RRS_REORDERING); + // { + // class ConstraintsReorderingHelper { + // void operator ()(dxQuickStepperStage4CallContext *stage4CallContext, unsigned int startIndex, unsigned int indicesCount) + // { + // IndexError *order = stage4CallContext->m_order + startIndex; + // + // for (unsigned int index = 1; index < indicesCount; ++index) { + // int swapIndex = dRandInt(index + 1); + // IndexError tmp = order[index]; + // order[index] = order[swapIndex]; + // order[swapIndex] = tmp; + // } + // } + // }; + + /* + * Full reorder needs to be done. + * Even though this contradicts the initial idea of moving dependent constraints + * to the order end the algorithm does not work the other way well. + * It looks like the iterative method needs a shake after it already found + * some initial approximations and those incurred errors help it to converge even better. + */ if (Atomics.ThrsafeExchange(stage4CallContext.m_SOR_reorderHeadTaken, 1) == 0) { // Process the head - dxQuickStepperLocalContext localContext = stage4CallContext.m_localContext; - ConstraintsReorderingHelper(stage4CallContext, 0, - localContext.m_m - localContext.m_valid_findices.get()); - } - if (Atomics.ThrsafeExchange(stage4CallContext.m_SOR_reorderTailTaken, 1) == 0) { - // Process the tail - dxQuickStepperLocalContext localContext = stage4CallContext.m_localContext; - ConstraintsReorderingHelper(stage4CallContext, localContext.m_m - localContext.m_valid_findices.get(), - localContext.m_valid_findices.get()); + final dxQuickStepperLocalContext localContext = stage4CallContext.m_localContext; + ConstraintsReorderingHelper(stage4CallContext, 0, localContext.m_m); } - } else { - // Revert to the normal order on the next step after the shuffling - dxQuickStepIsland_Stage4LCP_ReorderPrep(stage4CallContext); - } + } + // dIASSERT((RRS__MAX, true)); // A reference to RRS__MAX to be located by Find in Files } else { // Just return true and skip the randomization for the very first iteration } result = true; + // #else // #if CONSTRAINTS_REORDERING_METHOD != REORDERING_METHOD__BY_ERROR && CONSTRAINTS_REORDERING_METHOD != REORDERING_METHOD__RANDOMLY } else { - // TODO CHECK-TZ Confirm that this works, otherwise remove condition again - if (iteration == 0) { - result = true; // return true on 0th iteration to build dependency map for initial order - } + dIASSERT(iteration == 0); // The reordering request is only returned for the first iteration + result = true; } return result; } @@ -1966,6 +2055,23 @@ void dxQuickStepIsland_Stage4LCP_DependencyMapFromSavedLevelsReconstruction(dxQu AtomicInteger[] mi_levels = stage4CallContext.m_bi_links_or_mi_levels; AtomicInteger[] mi_links = stage4CallContext.m_mi_links; + // NOTE! + // OD: The mi_links array is not zero-filled before the reconstruction. + // Iteration ends with all the down links zeroed. And since down links + // are moved to the next level links when parent-child relations are established, + // the horizontal levels are properly terminated. + // The leaf nodes had their links zero-initialized initially + // and those zeros remain intact during the solving. This way the down links + // are properly terminated as well. + // This is very obscure and error prone and would need an assertion check at least + // but the simplest assertion approach I can imagine would be + // zero filling and building another tree with the memory buffer comparison afterwards. + // That would be stupid, obviously. + // + // NOTE! + // OD: This routine can be threaded. However having two threads messing + // in one integer array with random access and kicking each other memory lines + // out of cache would probably work worse than letting a single thread do the whole job. int m = localContext.m_m; for (int i = 0; i != m; ++i) { int currentLevelRoot = mi_levels[i].get(); @@ -1998,10 +2104,37 @@ public void run() { } dxQuickStepIsland_Stage4LCP_MTIteration(stage4CallContext, knownToBeCompletedLevel); } - - public static AtomicInteger mtIterations = new AtomicInteger(); - - private static + + public static AtomicInteger mtIterations = new AtomicInteger(); + + /* + * +0 +0 + * Root───┬─────────────────┬──... + * +1│ +1│ + * ┌┴┐+0 ┌─┐+0 . + * │A├─────┤B├─... + * └┬┘ └┬┘ + * +1│ +1│ + * ┌┴┐+0 . + * │C├─... + * └┬┘ + * +1│ + * . + * + * Lower tree levels depend on their parents. Same level nodes are independent with respect to each other. + * + * 1. B is linked in place of A + * 2. A is processed + * 3. C is inserted at the Root level + * + * The tree starts with a single child subtree at the root level ("down" link of slot #0 is used for that). + * Then, additional "C" nodes are added to the root level by building horizontal link via slots of + * their former parent "A"s that had become free. + * The "level" link of slot #0 is used to find the root level head. + * + * Since the tree is altered during iteration, mi_levels record each node parents so that the tree could be reconstructed. + */ + private static void dxQuickStepIsland_Stage4LCP_MTIteration(final dxQuickStepperStage4CallContext stage4CallContext, final int initiallyKnownToBeCompletedLevel) { mtIterations.incrementAndGet(); @@ -2247,6 +2380,8 @@ void dxQuickStepIsland_Stage4b(dxQuickStepperStage4CallContext stage4CallContext // they should not be used again. if (IsStage4bJointInfosIterationRequired(localContext)) { + DVector3 dataL = new DVector3(); //JVE__MAX; + DVector3 dataA = new DVector3(); //JVE__MAX; double[] Jcopy = localContext.m_Jcopy; double[] lambda = stage4CallContext.m_lambda; int[] mindex = localContext.m_mindex; @@ -2259,33 +2394,75 @@ void dxQuickStepIsland_Stage4b(dxQuickStepperStage4CallContext stage4CallContext int ji_step; while ((ji_step = Atomics.ThrsafeIncrementIntUpToLimit(stage4CallContext.m_ji_4b, nj_steps)) != nj_steps) { int ji = ji_step * step_size; - int lambdacurr = mindex[ji * 2]; - int Jcopycurr = mindex[ji * 2 + 1] * JCE__MAX; - int jicurr = ji; - int jiend = jicurr + Math.min(step_size, nj - ji); + final int jiend = ji + Math.min(step_size, nj - ji); - while (true) { - DxJoint joint = jointinfos[jicurr].joint; - int infom = jointinfos[jicurr].info.m; + //const dReal *Jcopycurr = Jcopy + (sizeint)mindex[ji].fbIndex * JCE__MAX; + int Jcopycurr = getFbIndex(mindex, ji) * JCE__MAX; + while (true) { // straightforward computation of joint constraint forces: // multiply related lambdas with respective J' block for joints // where feedback was requested - DJoint.DJointFeedback fb = joint.feedback; - if (fb != null) { - dAssertVec3Element(); // Assigning fb.f1/f2/t1/t2 inside multiply, see ODE code. - Multiply1_12q1 (fb.f1, fb.t1, Jcopy, Jcopycurr + JCE__J1_MIN, lambda, lambdacurr, infom); - if (joint.node[1].body != null) { - Multiply1_12q1 (fb.f2, fb.t2, Jcopy, Jcopycurr + JCE__J2_MIN, lambda, lambdacurr, infom); - } - Jcopycurr += infom * JCE__MAX; - } + final int fb_infom = getFbIndex(mindex, ji + 1) - getFbIndex(mindex, ji); + if (fb_infom != 0) { + dIASSERT(fb_infom == getMIndex(mindex, ji + 1) - getMIndex(mindex, ji)); + + //const dReal *lambdacurr = lambda + mindex[ji].mIndex; + final int lambdacurrOfs = getMIndex(mindex, ji); + DxJoint joint = jointinfos[ji].joint; + + // #ifdef WARM_STARTING + if (WARM_STARTING) { + memcpy(joint.lambda, 0, lambda, lambdacurrOfs, fb_infom); + } + // #endif + + DJoint.DJointFeedback fb = joint.feedback; + + dAssertVec3Element(); // ode4j specific assertion + if (joint.node[1].body != null) { + Multiply1_12q1 (dataL, dataA, Jcopy, Jcopycurr + JCE__J2_MIN, lambda, lambdacurrOfs, fb_infom); + dSASSERT(JCE__MAX == 12); + + // fb.f2[dSA_X] = data[JVE_LX]; + // fb->f2[dSA_Y] = data[JVE_LY]; + // fb->f2[dSA_Z] = data[JVE_LZ]; + // fb->t2[dSA_X] = data[JVE_AX]; + // fb->t2[dSA_Y] = data[JVE_AY]; + // fb->t2[dSA_Z] = data[JVE_AZ]; + fb.f2.set(dataL); + fb.t2.set(dataA); + } + + Multiply1_12q1 (dataL, dataA, Jcopy, Jcopycurr + JCE__J1_MIN, lambda, lambdacurrOfs, fb_infom); + dSASSERT(JCE__MAX == 12); + + // fb->f1[dSA_X] = data[JVE_LX]; + // fb->f1[dSA_Y] = data[JVE_LY]; + // fb->f1[dSA_Z] = data[JVE_LZ]; + // fb->t1[dSA_X] = data[JVE_AX]; + // fb->t1[dSA_Y] = data[JVE_AY]; + // fb->t1[dSA_Z] = data[JVE_AZ]; + fb.f2.set(dataL); + fb.t2.set(dataA); + + Jcopycurr += fb_infom * JCE__MAX; + } + else { + // #ifdef WARM_STARTING + if (WARM_STARTING) { + int lambdacurrOfs = getMIndex(mindex, ji); + int infom = getMIndex(mindex, ji + 1) - getMIndex(mindex, ji); + DxJoint joint = jointinfos[ji].joint; + memcpy(joint.lambda, 0, lambda, lambdacurrOfs, infom); + // #endif + } + } - if (++jicurr == jiend) { - break; - } - lambdacurr += infom; - } + if (++ji == jiend) { + break; + } + } } } } diff --git a/core/src/main/java/org/ode4j/ode/internal/DxSAPSpace.java b/core/src/main/java/org/ode4j/ode/internal/DxSAPSpace.java index 9e92e879..f5c48aa7 100644 --- a/core/src/main/java/org/ode4j/ode/internal/DxSAPSpace.java +++ b/core/src/main/java/org/ode4j/ode/internal/DxSAPSpace.java @@ -261,7 +261,13 @@ void remove( DxGeom g ) GEOM_SET_GEOM_IDX(g,GEOM_INVALID_IDX); GeomList.remove( geomSize-1 ); } - + + //g->tome_ex = 0; + g._qtIdxEx = null; + // dUASSERT((g->next_ex = 0, true), "Needed for an assertion check only"); + g.setNextEx(null); + dUASSERT(true, "Needed for an assertion check only"); + super.remove(g); } diff --git a/core/src/main/java/org/ode4j/ode/internal/DxWorld.java b/core/src/main/java/org/ode4j/ode/internal/DxWorld.java index 62ca30f3..0364870a 100644 --- a/core/src/main/java/org/ode4j/ode/internal/DxWorld.java +++ b/core/src/main/java/org/ode4j/ode/internal/DxWorld.java @@ -281,50 +281,50 @@ int dWorldGetStepIslandsProcessingMaxThreadCount() return islands_max_threads; } - /** - * - * @param policyinfo - * @return x - * @deprecated Not used anymore: remove TODO - */ - @Deprecated - boolean dWorldSetStepMemoryReservationPolicy(final DWorldStepReserveInfo policyinfo) - { - dUASSERT (policyinfo==null || (policyinfo.struct_size >= DxUtil.sizeof(policyinfo) && policyinfo.reserve_factor >= 1.0f), "Bad policy info"); - - boolean result = false; - - //(TZ) DxStepWorkingMemory wmem = policyinfo!=null ? AllocateOnDemand(this.wmem) : this.wmem; - DxStepWorkingMemory wmem; - if (policyinfo!=null) { - if (this.wmem == null) { - this.wmem = new DxStepWorkingMemory(); - } - wmem = this.wmem; - } else { - wmem = this.wmem; - } - - if (wmem!=null) - { - if (policyinfo!=null) - { - wmem.SetMemoryReserveInfo(policyinfo.reserve_factor, policyinfo.reserve_minimum); - result = wmem.GetMemoryReserveInfo() != null; - } - else - { - wmem.ResetMemoryReserveInfoToDefault(); - result = true; - } - } - else if (policyinfo==null) - { - result = true; - } - - return result; - } + // /** + // * + // * @param policyinfo + // * @return x + // * @deprecated Not used anymore + // */ + // @Deprecated + // boolean dWorldSetStepMemoryReservationPolicy(final DWorldStepReserveInfo policyinfo) + // { + // dUASSERT (policyinfo==null || (policyinfo.struct_size >= DxUtil.sizeof(policyinfo) && policyinfo.reserve_factor >= 1.0f), "Bad policy info"); + // + // boolean result = false; + // + // //(TZ) DxStepWorkingMemory wmem = policyinfo!=null ? AllocateOnDemand(this.wmem) : this.wmem; + // DxStepWorkingMemory wmem; + // if (policyinfo!=null) { + // if (this.wmem == null) { + // this.wmem = new DxStepWorkingMemory(); + // } + // wmem = this.wmem; + // } else { + // wmem = this.wmem; + // } + // + // if (wmem!=null) + // { + // if (policyinfo!=null) + // { + // wmem.SetMemoryReserveInfo(policyinfo.reserve_factor, policyinfo.reserve_minimum); + // result = wmem.GetMemoryReserveInfo() != null; + // } + // else + // { + // wmem.ResetMemoryReserveInfoToDefault(); + // result = true; + // } + // } + // else if (policyinfo==null) + // { + // result = true; + // } + // + // return result; + // } boolean dWorldStep (double stepsize) diff --git a/core/src/main/java/org/ode4j/ode/internal/ErrorHdl.java b/core/src/main/java/org/ode4j/ode/internal/ErrorHdl.java index 376591b7..9a1ed091 100644 --- a/core/src/main/java/org/ode4j/ode/internal/ErrorHdl.java +++ b/core/src/main/java/org/ode4j/ode/internal/ErrorHdl.java @@ -119,10 +119,11 @@ public static void dError (int num, String msg, Object ... ap) { // va_start (ap,msg); if (error_function != null) { error_function.call (num,msg,ap); - } - logger.error("ODE Error " + num + ": " + msg, ap); + } else { + logger.error("ODE Error " + num + ": " + msg, ap); + } throw new RuntimeException("#"+num + ": " + msg); - //System.exit (1); + //exit (1); } @@ -137,8 +138,9 @@ public static void dDebug (int num, String msg, Object ... ap) { // va_start (ap,msg); if (debug_function != null) { debug_function.call (num,msg,ap); + } else { + logger.debug("ODE INTERNAL ERROR " + " " + num + ": " + msg, ap); } - logger.debug("ODE INTERNAL ERROR " + " " + num + ": " + msg, ap); // *((char *)0) = 0; ... commit SEGVicide //abort(); throw new RuntimeException("#"+num + ": " + String.format(msg, ap)); diff --git a/core/src/main/java/org/ode4j/ode/internal/FastDot.java b/core/src/main/java/org/ode4j/ode/internal/FastDot.java index b5adb40a..e7b4e643 100644 --- a/core/src/main/java/org/ode4j/ode/internal/FastDot.java +++ b/core/src/main/java/org/ode4j/ode/internal/FastDot.java @@ -103,7 +103,8 @@ public static double dDot (final double[] a, int aOfs, } //template - public static double dxtDot (final double[] a, int a_pos, final double[] b, int b_pos, int n, int b_stride) + public static double calculateLargeVectorDot ( + final double[] a, int a_pos, final double[] b, int b_pos, int n, final int b_stride) { double sum = 0; int a_end = a_pos + (n & (int)(~3)); @@ -122,5 +123,4 @@ public static double dxtDot (final double[] a, int a_pos, final double[] b, int } return sum; } - } \ No newline at end of file diff --git a/core/src/main/java/org/ode4j/ode/internal/FastLDLT.java b/core/src/main/java/org/ode4j/ode/internal/FastLDLTFactor.java similarity index 94% rename from core/src/main/java/org/ode4j/ode/internal/FastLDLT.java rename to core/src/main/java/org/ode4j/ode/internal/FastLDLTFactor.java index 4ed30bb1..1da4f375 100644 --- a/core/src/main/java/org/ode4j/ode/internal/FastLDLT.java +++ b/core/src/main/java/org/ode4j/ode/internal/FastLDLTFactor.java @@ -30,7 +30,7 @@ /** * Implementation of FastDLTImpl */ -public class FastLDLT { +public class FastLDLTFactor { //#ifndef _ODE_FASTLDLT_IMPL_H_ //#define _ODE_FASTLDLT_IMPL_H_ @@ -48,7 +48,7 @@ public class FastLDLT { // template - public static void dxtFactorLDLT(double[] A, double[] d, int rowCount, int rowSkip, int d_stride) { + public static void factorMatrixAsLDLT(double[] A, double[] d, int rowCount, int rowSkip, int d_stride) { if (rowCount < 1) return; final int lastRowIndex = rowCount - 1; @@ -60,10 +60,10 @@ public static void dxtFactorLDLT(double[] A, double[] d, int rowCount, int rowSk for (; blockStartRow < lastRowIndex; subsequentPass = true, ARow_pos += 2 * rowSkip, blockStartRow += 2) { if (subsequentPass) { /* solve L*(D*l)=a, l is scaled elements in 2 x i block at A(i,0) */ - dxSolveL1_2(A, A, ARow_pos, blockStartRow, rowSkip); - dxScaleAndFactorizeL1_2(A, ARow_pos, d, blockStartRow, rowSkip, d_stride); + solveL1Stripe_2(A, A, ARow_pos, blockStartRow, rowSkip); + scaleAndFactorizeL1Stripe_2(A, ARow_pos, d, blockStartRow, rowSkip, d_stride); } else { - dxScaleAndFactorizeFirstL1Row_2(A, ARow_pos, d, rowSkip, d_stride); + scaleAndFactorizeL1FirstRowStripe_2(A, ARow_pos, d, rowSkip, d_stride); } /* done factorizing 2 x 2 block */ } @@ -71,10 +71,10 @@ public static void dxtFactorLDLT(double[] A, double[] d, int rowCount, int rowSk /* compute the (less than 2) rows at the bottom */ if (!subsequentPass || blockStartRow == lastRowIndex) { if (subsequentPass) { - dxSolveL1_1(A, A, ARow_pos, blockStartRow, rowSkip); - dxScaleAndFactorizeL1_1(A, ARow_pos, d, blockStartRow, d_stride); + solveStripeL1_1(A, A, ARow_pos, blockStartRow, rowSkip); + scaleAndFactorizeL1Stripe_1(A, ARow_pos, d, blockStartRow, d_stride); } else { - dxScaleAndFactorizeFirstL1Row_1(A, ARow_pos, d, d_stride); + scaleAndFactorizeL1FirstRowStripe_1(A, ARow_pos, d, d_stride); } /* done factorizing 1 x 1 block */ } @@ -89,7 +89,7 @@ public static void dxtFactorLDLT(double[] A, double[] d, int rowCount, int rowSk * this processes blocks of 2*2. * if this is in the factorizer source file, n must be a multiple of 2. */ - static void dxSolveL1_2(final double[] L, double[] B, int bPos, int rowCount, int rowSkip) { + static void solveL1Stripe_2(final double[] L, double[] B, int bPos, int rowCount, int rowSkip) { dIASSERT(rowCount != 0); dIASSERT(rowCount % 2 == 0); @@ -229,7 +229,7 @@ static void dxSolveL1_2(final double[] L, double[] B, int bPos, int rowCount, in } //template - private static void dxScaleAndFactorizeL1_2(double[] ARow, int aPos, double[] d, int factorizationRow, int rowSkip, int d_stride) { + private static void scaleAndFactorizeL1Stripe_2(double[] ARow, int aPos, double[] d, int factorizationRow, int rowSkip, int d_stride) { dIASSERT(factorizationRow != 0); dIASSERT(factorizationRow % 2 == 0); @@ -346,7 +346,7 @@ private static void dxScaleAndFactorizeL1_2(double[] ARow, int aPos, double[] d, } //template - private static void dxScaleAndFactorizeFirstL1Row_2(double[] ARow, int aPos, double[] d, int rowSkip, int d_stride) { + private static void scaleAndFactorizeL1FirstRowStripe_2(double[] ARow, int aPos, double[] d, int rowSkip, int d_stride) { int ptrAElement = aPos; //ARow; double[] ptrDElement = d; @@ -380,7 +380,7 @@ private static void dxScaleAndFactorizeFirstL1Row_2(double[] ARow, int aPos, dou * this processes blocks of 2*2. * if this is in the factorizer source file, n must be a multiple of 2. */ - static void dxSolveL1_1(final double[] L, double[] B, int bPos, int rowCount, int rowSkip) { + static void solveStripeL1_1(final double[] L, double[] B, int bPos, int rowCount, int rowSkip) { dIASSERT(rowCount != 0); dIASSERT(rowCount % 2 == 0); @@ -491,7 +491,7 @@ static void dxSolveL1_1(final double[] L, double[] B, int bPos, int rowCount, in } //template - private static void dxScaleAndFactorizeL1_1(double[] ARow, int APos, double[] d, int factorizationRow, int d_stride) { + private static void scaleAndFactorizeL1Stripe_1(double[] ARow, int APos, double[] d, int factorizationRow, int d_stride) { int ptrAElement = APos;//ARow; int ptrDElement = 0;//d; @@ -561,7 +561,7 @@ private static void dxScaleAndFactorizeL1_1(double[] ARow, int APos, double[] d, } //template - private static void dxScaleAndFactorizeFirstL1Row_1(double[] ARow, int APos, double[] d, int d_stride) { + private static void scaleAndFactorizeL1FirstRowStripe_1(double[] ARow, int APos, double[] d, int d_stride) { int ptrAElement = APos; //ARow; double[] ptrDElement = d; diff --git a/core/src/main/java/org/ode4j/ode/internal/FastLDLTSolve.java b/core/src/main/java/org/ode4j/ode/internal/FastLDLTSolve.java new file mode 100644 index 00000000..5159a8e5 --- /dev/null +++ b/core/src/main/java/org/ode4j/ode/internal/FastLDLTSolve.java @@ -0,0 +1,53 @@ +/************************************************************************* + * * + * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith. * + * All rights reserved. Email: russ@q12.org Web: www.q12.org * + * Open Dynamics Engine 4J, Copyright (C) 2009-2021 Tilmann Zaeschke * + * All rights reserved. Email: ode4j@gmx.de Web: www.ode4j.org * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of EITHER: * + * (1) The GNU Lesser General Public License as published by the Free * + * Software Foundation; either version 2.1 of the License, or (at * + * your option) any later version. The text of the GNU Lesser * + * General Public License is included with this library in the * + * file LICENSE.TXT. * + * (2) The BSD-style license that is included with this library in * + * the file LICENSE-BSD.TXT. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files * + * LICENSE.TXT and LICENSE-BSD.TXT for more details. * + * * + *************************************************************************/ +package org.ode4j.ode.internal; + + +import static org.ode4j.ode.internal.Common.dAASSERT; +import static org.ode4j.ode.internal.FastLSolve.solveL1Straight; +import static org.ode4j.ode.internal.FastLTSolve.solveL1Transposed; +import static org.ode4j.ode.internal.FastVecScale.scaleLargeVector; + +public class FastLDLTSolve { + + + // #ifndef _ODE_MATRIX_IMPL_H_ + // #define _ODE_MATRIX_IMPL_H_ + + + // template + // void solveEquationSystemWithLDLT (const dReal *L, const dReal *d, dReal *b, unsigned rowCount, unsigned rowSkip) + public static void solveEquationSystemWithLDLT(final double[] L, double[] dArray, int dPos, double[] bArray, int bPos, int rowCount, int rowSkip, int d_stride, int b_stride) { + dAASSERT(L != null); + dAASSERT(dArray != null); + dAASSERT(bArray != null); + dAASSERT(rowCount > 0); + dAASSERT(rowSkip >= rowCount); + + solveL1Straight(L, bArray, bPos, rowCount, rowSkip, b_stride); + scaleLargeVector(bArray, bPos, dArray, dPos, rowCount, b_stride, d_stride); + solveL1Transposed(L, bArray, bPos, rowCount, rowSkip, b_stride); + } + +} diff --git a/core/src/main/java/org/ode4j/ode/internal/FastLSolve.java b/core/src/main/java/org/ode4j/ode/internal/FastLSolve.java index 6bd8f079..75f039f5 100644 --- a/core/src/main/java/org/ode4j/ode/internal/FastLSolve.java +++ b/core/src/main/java/org/ode4j/ode/internal/FastLSolve.java @@ -43,9 +43,9 @@ public class FastLSolve { * if this is in the factorizer source file, n must be a multiple of 4. */ //template - public static void dxtSolveL1 (final double[] L, double[] B, final int BPos, int rowCount, int rowSkip, int b_stride) + public static void solveL1Straight (final double[] L, double[] B, final int BPos, int rowCount, int rowSkip, int b_stride) { - int LPos = 0; + final int LPos = 0; dIASSERT(rowCount != 0); /* compute all 4 x 1 blocks of X */ diff --git a/core/src/main/java/org/ode4j/ode/internal/FastLTSolve.java b/core/src/main/java/org/ode4j/ode/internal/FastLTSolve.java index cb3c66f5..6afed4fc 100644 --- a/core/src/main/java/org/ode4j/ode/internal/FastLTSolve.java +++ b/core/src/main/java/org/ode4j/ode/internal/FastLTSolve.java @@ -40,7 +40,7 @@ public class FastLTSolve { * this processes blocks of 4. */ //template - public static void dxtSolveL1T(final double[] L, double[] B, int BPos, int rowCount, int rowSkip, int b_stride) + public static void solveL1Transposed(final double[] L, double[] B, int BPos, int rowCount, int rowSkip, int b_stride) { int LPos = 0; dIASSERT(rowCount != 0); diff --git a/core/src/main/java/org/ode4j/ode/internal/FastVecScale.java b/core/src/main/java/org/ode4j/ode/internal/FastVecScale.java new file mode 100644 index 00000000..d13ea1b6 --- /dev/null +++ b/core/src/main/java/org/ode4j/ode/internal/FastVecScale.java @@ -0,0 +1,89 @@ +/************************************************************************* + * * + * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith. * + * All rights reserved. Email: russ@q12.org Web: www.q12.org * + * Open Dynamics Engine 4J, Copyright (C) 2009-2021 Tilmann Zaeschke * + * All rights reserved. Email: ode4j@gmx.de Web: www.ode4j.org * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of EITHER: * + * (1) The GNU Lesser General Public License as published by the Free * + * Software Foundation; either version 2.1 of the License, or (at * + * your option) any later version. The text of the GNU Lesser * + * General Public License is included with this library in the * + * file LICENSE.TXT. * + * (2) The BSD-style license that is included with this library in * + * the file LICENSE-BSD.TXT. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files * + * LICENSE.TXT and LICENSE-BSD.TXT for more details. * + * * + *************************************************************************/ +package org.ode4j.ode.internal; + + +import static org.ode4j.ode.internal.Common.dAASSERT; + +public class FastVecScale { + + // #ifndef _ODE_FASTVECSCALE_IMPL_H_ + // #define _ODE_FASTVECSCALE_IMPL_H_ + + + + // matrix_impl.h + + // template + // void scaleLargeVector (dReal *aStart, const dReal *dStart, unsigned elementCount) + public static void scaleLargeVector(double[] A, final int aStart, final double[] D, int dStart, int elementCount, int a_stride, int d_stride) { + dAASSERT(aStart >= 0 && dStart >= 0 && elementCount >= 0); + + final int step = 4; + // dReal *ptrA = aStart; + // const dReal *ptrD = dStart; + // const dReal *const dStepsEnd = dStart + (size_t)(elementCount & ~(step - 1)) * d_stride; + int ptrA = aStart; + int ptrD = dStart; + final int dStepsEnd = dStart + (elementCount & ~(step - 1)) * d_stride; + + for (; ptrD != dStepsEnd; ptrA += step * a_stride, ptrD += step * d_stride) { + double a0 = A[ptrA + 0], a1 = A[ptrA + 1 * a_stride], a2 = A[ptrA + 2 * a_stride], a3 = A[ptrA + 3 * a_stride]; + double d0 = D[ptrD + 0], d1 = D[ptrD + 1 * d_stride], d2 = D[ptrD + 2 * d_stride], d3 = D[ptrD + 3 * d_stride]; + a0 *= d0; + a1 *= d1; + a2 *= d2; + a3 *= d3; + A[ptrA + 0] = a0; + A[ptrA + 1 * a_stride] = a1; + A[ptrA + 2 * a_stride] = a2; + A[ptrA + 3 * a_stride] = a3; + //dSASSERT(step == 4); + } + + switch (elementCount & (step - 1)) { + case 3: { + double a2 = A[ptrA + 2 * a_stride]; + double d2 = D[ptrD + 2 * d_stride]; + A[ptrA + 2 * a_stride] = a2 * d2; + // break; -- proceed to case 2 + } + + case 2: { + double a1 = A[ptrA + 1 * a_stride]; + double d1 = D[ptrD + 1 * d_stride]; + A[ptrA + 1 * a_stride] = a1 * d1; + // break; -- proceed to case 1 + } + + case 1: { + double a0 = A[ptrA + 0]; + double d0 = D[ptrD + 0]; + A[ptrA + 0] = a0 * d0; + break; + } + } + } + +} diff --git a/core/src/main/java/org/ode4j/ode/internal/Matrix.java b/core/src/main/java/org/ode4j/ode/internal/Matrix.java index 814f8036..c9aaafcf 100644 --- a/core/src/main/java/org/ode4j/ode/internal/Matrix.java +++ b/core/src/main/java/org/ode4j/ode/internal/Matrix.java @@ -24,8 +24,9 @@ *************************************************************************/ package org.ode4j.ode.internal; -import static org.ode4j.ode.internal.FastLSolve.dxtSolveL1; -import static org.ode4j.ode.internal.FastLTSolve.dxtSolveL1T; +import static org.ode4j.ode.internal.FastLDLTSolve.solveEquationSystemWithLDLT; +import static org.ode4j.ode.internal.FastLSolve.solveL1Straight; +import static org.ode4j.ode.internal.FastLTSolve.solveL1Transposed; import static org.ode4j.ode.internal.cpp4j.Cstring.memcpy; import static org.ode4j.ode.internal.cpp4j.Cstring.memmove; @@ -348,8 +349,7 @@ public static void dMultiply1(double[] A, final double[] B, * @param B B * @param C C */ - public static void dMultiply2(DMatrix3 A, final DMatrix3C B, - final DMatrix3C C) { + public static void dMultiply2(DMatrix3 A, final DMatrix3C B, final DMatrix3C C) { //dMultiply2(A.v, ((DMatrix3) B).v, ((DMatrix3) C).v, 3, 3, 3); dMultiply0(A, B, C.reTranspose()); } @@ -368,8 +368,7 @@ public static void dMultiply2(DMatrix3 A, final DMatrix3C B, * @param B B * @param C C */ - public static void dMultiply2(DVector3 A, final DMatrix3C B, - final DVector3C C) { + public static void dMultiply2(DVector3 A, final DMatrix3C B, final DVector3C C) { //dMultiply2(A.v, ((DMatrix3) B).v, ((DVector3) C).v, 3, 3, 1); //TZ: this is equal to dMultiply0(...) !!! A.set0( B.get00()*C.get0() + B.get01()*C.get1() + B.get02()*C.get2() ); @@ -393,8 +392,7 @@ public static void dMultiply2(DVector3 A, final DMatrix3C B, * @param r r * */ - public static void dMultiply2(double[] A, final double[] B, - final double[] C, int p, int q, int r) { + public static void dMultiply2(double[] A, final double[] B, final double[] C, int p, int q, int r) { dAASSERT(p > 0 && q > 0 && r > 0); final int rskip = dPAD(r); final int qskip = dPAD(q); @@ -412,37 +410,6 @@ public static void dMultiply2(double[] A, final double[] B, A[a] = sum;//(*a) = sum; } } - //TODO remove is from 0.11.1 -// int i, j, k, z, rpad, qskip; -// double sum; -// // final double[] bb,cc; -// // TZ: -// int aPos = 0, bPos, cPos; -// // dAASSERT (A, B , C); -// dAASSERT(p > 0 && q > 0 && r > 0); -// rpad = dPAD(r) - r; -// qskip = dPAD(q); -// // bb = B; -// bPos = 0; -// for (i = p; i > 0; i--) { -// // cc = C; -// cPos = 0; -// for (j = r; j > 0; j--) { -// z = 0; -// sum = 0; -// // for (k=q; k>0; k--,z++) sum += bb[z] * cc[z]; -// for (k = q; k > 0; k--, z++) -// sum += B[bPos + z] * C[cPos + z]; -// // *(A++) = sum; -// A[aPos++] = sum; -// // cc += qskip; -// cPos += qskip; -// } -// // A += rpad; -// aPos += rpad; -// // bb += qskip; -// bPos += qskip; -// } } /** @@ -713,18 +680,6 @@ public static boolean dIsPositiveDefinite(DMatrix3C A) { return dFactorCholesky(A.clone()); } -// /** -// * this has been replaced by a faster version void dSolveL1T (const double -// * *L, double *b, int n, int nskip) { int i,j; dAASSERT (L && b && n >= 0 && -// * nskip >= n); double sum; for (i=n-2; i>=0; i--) { sum = 0; for (j=i+1; -// * j - // void dxtVectorScale (dReal *aStart, const dReal *dStart, unsigned elementCount) - public static void dxtVectorScale(double[] A, final int aStart, final double[] D, int dStart, int elementCount, int a_stride, int d_stride) { - dAASSERT(aStart >= 0 && dStart >= 0 && elementCount >= 0); - - final int step = 4; - // dReal *ptrA = aStart; - // const dReal *ptrD = dStart; - // const dReal *const dStepsEnd = dStart + (size_t)(elementCount & ~(step - 1)) * d_stride; - int ptrA = aStart; - int ptrD = dStart; - final int dStepsEnd = dStart + (elementCount & ~(step - 1)) * d_stride; - - for (; ptrD != dStepsEnd; ptrA += step * a_stride, ptrD += step * d_stride) { - double a0 = A[ptrA + 0], a1 = A[ptrA + 1 * a_stride], a2 = A[ptrA + 2 * a_stride], a3 = A[ptrA + 3 * a_stride]; - double d0 = D[ptrD + 0], d1 = D[ptrD + 1 * d_stride], d2 = D[ptrD + 2 * d_stride], d3 = D[ptrD + 3 * d_stride]; - a0 *= d0; - a1 *= d1; - a2 *= d2; - a3 *= d3; - A[ptrA + 0] = a0; - A[ptrA + 1 * a_stride] = a1; - A[ptrA + 2 * a_stride] = a2; - A[ptrA + 3 * a_stride] = a3; - //dSASSERT(step == 4); - } - - switch (elementCount & (step - 1)) { - case 3: { - double a2 = A[ptrA + 2 * a_stride]; - double d2 = D[ptrD + 2 * d_stride]; - A[ptrA + 2 * a_stride] = a2 * d2; - // break; -- proceed to case 2 - } - - case 2: { - double a1 = A[ptrA + 1 * a_stride]; - double d1 = D[ptrD + 1 * d_stride]; - A[ptrA + 1 * a_stride] = a1 * d1; - // break; -- proceed to case 1 - } - - case 1: { - double a0 = A[ptrA + 0]; - double d0 = D[ptrD + 0]; - A[ptrA + 0] = a0 * d0; - break; - } - } - } - - - // template - // void dxtSolveLDLT (const dReal *L, const dReal *d, dReal *b, unsigned rowCount, unsigned rowSkip) - public static void dxtSolveLDLT(final double[] L, double[] dArray, int dPos, double[] bArray, int bPos, int rowCount, int rowSkip, int d_stride, int b_stride) { - dAASSERT(L != null); - dAASSERT(dArray != null); - dAASSERT(bArray != null); - dAASSERT(rowCount > 0); - dAASSERT(rowSkip >= rowCount); - - dxtSolveL1(L, bArray, bPos, rowCount, rowSkip, b_stride); - dxtVectorScale(bArray, bPos, dArray, dPos, rowCount, b_stride, d_stride); - dxtSolveL1T(L, bArray, bPos, rowCount, rowSkip, b_stride); - } - public static double dxDot(final double[] a, final double[] b, int n) { return dDot(a, b, n); } diff --git a/core/src/main/java/org/ode4j/ode/internal/Step.java b/core/src/main/java/org/ode4j/ode/internal/Step.java index 9eea11a6..062ce09c 100644 --- a/core/src/main/java/org/ode4j/ode/internal/Step.java +++ b/core/src/main/java/org/ode4j/ode/internal/Step.java @@ -1582,23 +1582,14 @@ void dxStepIsland_Stage3(dxStepperStage3CallContext stage3CallContext) stage3CallContext = null; // WARNING! stage3CallContext is not valid after this point! dIVERIFY(stage3CallContext == null); // To suppress unused variable assignment warnings - double[] invI = localContext.m_invI; - dJointWithInfo1[] jointinfosA = localContext.m_jointinfosA; - int jointinfosP = localContext.m_jointinfosOfs; - int nj = localContext.m_nj; int m = localContext.m_m; int nub = localContext.m_nub; //const unsigned int *mindex = localContext.m_mindex; int[] findex = localContext.m_findex; - double[] J = localContext.m_J; double[] A = localContext.m_A; double[] pairsRhsLambda = localContext.m_pairsRhsCfm; // Reuse cfm buffer for lambdas as the former values are not needed any more double[] pairsLoHi = localContext.m_pairsLoHi; - DxBody[] bodyA = callContext.m_islandBodiesStartA(); - int bodyP = callContext.m_islandBodiesStartOfs(); - int nb = callContext.m_islandBodiesCount(); - if (m > 0) { BlockPointer lcpstate = memarena.BEGIN_STATE_SAVE(); { diff --git a/core/src/main/java/org/ode4j/ode/internal/joints/DxJoint.java b/core/src/main/java/org/ode4j/ode/internal/joints/DxJoint.java index 1a33f880..5c623322 100644 --- a/core/src/main/java/org/ode4j/ode/internal/joints/DxJoint.java +++ b/core/src/main/java/org/ode4j/ode/internal/joints/DxJoint.java @@ -409,6 +409,13 @@ void setBall2(DxJoint joint, double fps, double erp, int rowskip, double[] J1A, dCalcVectorCross3(J1A, J1Ofs + rowskip + JointEnums.GI2__JA_MIN, a1, q1); dCalcVectorCross3(J1A, J1Ofs + 2 * rowskip + JointEnums.GI2__JA_MIN, a1, q2); + DxBody b0 = node[0].body; + dAddVectors3(a1, a1, b0.posr().pos()); + + // set right hand side - measure error along (axis,q1,q2) + double k1 = fps * erp1; + double k = fps * erp; + DxBody b1 = joint.node[1].body; if (b1 != null) { dCopyNegatedVector3(J2A, J2Ofs + JointEnums.GI2__JL_MIN, axis); @@ -422,16 +429,7 @@ void setBall2(DxJoint joint, double fps, double erp, int rowskip, double[] J1A, dCalcVectorCross3(J2A, J2Ofs + 2 * rowskip + JointEnums.GI2__JA_MIN, q2, a2); //== dCalcVectorCross3( J2 + // 2 * rowskip + JointEnums.GI2__J2A_MIN, a2, q2 ); dNegateVector3( J2 + 2 * rowskip + JointEnums // .GI2__J2A_MIN ); - } - // set right hand side - measure error along (axis,q1,q2) - double k1 = fps * erp1; - double k = fps * erp; - - DxBody b0 = joint.node[0].body; - dAddVectors3(a1, a1, b0.posr().pos()); - - if (b1 != null) { dAddVectors3(a2, a2, b1.posr().pos()); DVector3 a2_minus_a1 = new DVector3(); @@ -788,7 +786,7 @@ public void dJointAttach (DBody b1, DBody b2) // check if the joint can not be attached to just one body dUASSERT (!((flags & dJOINT_TWOBODIES) != 0 && - ((body1 != null) ^ (body2 != null))), + ((body1 != null) != (body2 != null))), "joint can not be attached to just one body"); // remove any existing body attachments diff --git a/core/src/main/java/org/ode4j/ode/internal/joints/DxJointContact.java b/core/src/main/java/org/ode4j/ode/internal/joints/DxJointContact.java index 57dfa4ba..774e56dc 100644 --- a/core/src/main/java/org/ode4j/ode/internal/joints/DxJointContact.java +++ b/core/src/main/java/org/ode4j/ode/internal/joints/DxJointContact.java @@ -80,40 +80,69 @@ void getSureMaxInfo( SureMaxInfo info ) // make sure mu's >= 0, then calculate number of constraint rows and number // of unbounded rows. int m = 1, nub = 0; - boolean roll = (contact.surface.mode&OdeConstants.dContactRolling)!=0; - - if ( contact.surface.mu < 0 ) contact.surface.mu = 0; - - // Anisotropic sliding and rolling and spinning friction - if ( (contact.surface.mode & OdeConstants.dContactAxisDep) != 0 ) - { - if ( contact.surface.mu2 < 0 ) contact.surface.mu2 = 0; - if ( contact.surface.mu > 0 ) m++; - if ( contact.surface.mu2 > 0 ) m++; - if ( contact.surface.mu == dInfinity ) nub ++; - if ( contact.surface.mu2 == dInfinity ) nub ++; - if (roll) { - if ( contact.surface.rho < 0 ) contact.surface.rho = 0; - else m++; - if ( contact.surface.rho2 < 0 ) contact.surface.rho2 = 0; - else m++; - if ( contact.surface.rhoN < 0 ) contact.surface.rhoN = 0; - else m++; - - if ( contact.surface.rho == dInfinity ) nub++; - if ( contact.surface.rho2 == dInfinity ) nub++; - if ( contact.surface.rhoN == dInfinity ) nub++; - } + + // Anisotropic sliding and rolling and spinning friction + if ((contact.surface.mode & dContactAxisDep) != 0) { + if (contact.surface.mu < 0) { + contact.surface.mu = 0; + } + else if (contact.surface.mu > 0) { + if (contact.surface.mu == dInfinity) { nub++; } + m++; + } + + if (contact.surface.mu2 < 0) { + contact.surface.mu2 = 0; + } + else if (contact.surface.mu2 > 0) { + if (contact.surface.mu2 == dInfinity) { nub++; } + m++; + } + + if ((contact.surface.mode & dContactRolling) != 0) { + if (contact.surface.rho < 0) { + contact.surface.rho = 0; + } + else { + if (contact.surface.rho == dInfinity) { nub++; } + m++; + } + + if (contact.surface.rho2 < 0) { + contact.surface.rho2 = 0; + } + else { + if (contact.surface.rho2 == dInfinity) { nub++; } + m++; + } + + if (contact.surface.rhoN < 0) { + contact.surface.rhoN = 0; + } + else { + if (contact.surface.rhoN == dInfinity) { nub++; } + m++; + } + } } - else - { - if ( contact.surface.mu > 0 ) m += 2; - if ( contact.surface.mu == dInfinity ) nub += 2; - if (roll) { - if ( contact.surface.rho < 0 ) contact.surface.rho = 0; - else m+=3; - if ( contact.surface.rho == dInfinity ) nub += 3; - } + else { + if (contact.surface.mu < 0) { + contact.surface.mu = 0; + } + else if (contact.surface.mu > 0) { + if (contact.surface.mu == dInfinity) { nub += 2; } + m += 2; + } + + if ((contact.surface.mode & dContactRolling) != 0) { + if (contact.surface.rho < 0) { + contact.surface.rho = 0; + } + else { + if (contact.surface.rho == dInfinity) { nub += 3; } + m += 3; + } + } } the_m = m; @@ -131,33 +160,8 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A final int ROW_NORMAL = 0; final int ROW__OPTIONAL_MIN = 1; - // get normal, with sign adjusted for body1/body2 polarity - DVector3 normal = new DVector3(); - if (isFlagsReverse()) { - dCopyNegatedVector3(normal, contact.geom.normal); - } else { - dCopyVector3(normal, contact.geom.normal); - } - - // c1,c2 = contact points with respect to body PORs - DVector3 c1 = new DVector3(), c2 = new DVector3(0, 0, 0); - - DxBody b0 = node[0].body; - dSubtractVectors3(c1, contact.geom.pos, b0.posr().pos()); - - // set Jacobian for normal - dCopyVector3(J1A, J1Ofs + ROW_NORMAL * rowskip + GI2__JL_MIN, normal); - dCalcVectorCross3(J1A, J1Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, c1, normal); - - DxBody b1 = node[1].body; - if (b1 != null) { - dSubtractVectors3(c2, contact.geom.pos, b1.posr().pos()); - dCopyNegatedVector3(J2A, J2Ofs + ROW_NORMAL * rowskip + GI2__JL_MIN, normal); - dCalcVectorCross3(J2A, J2Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, normal, c2); //== dCalcVectorCross3( - // J2A, J2Ofs + GI2__JA_MIN, c2, normal ); dNegateVector3( J2A, J2Ofs + GI2__JA_MIN ); - } - final int surface_mode = contact.surface.mode; + // set right hand side and cfm value for normal double erp = (surface_mode & dContactSoftERP) != 0 ? contact.surface.soft_erp : worldERP; double k = worldFPS * erp; @@ -167,28 +171,59 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A double motionN = (surface_mode & dContactMotionN) != 0 ? contact.surface.motionN : 0.0; final double pushout = k * depth + motionN; + + boolean apply_bounce = (surface_mode & dContactBounce) != 0 && contact.surface.bounce_vel >= 0; + double outgoing = 0; + // note: this cap should not limit bounce velocity final double maxvel = world.contactp.max_vel; double c = pushout > maxvel ? maxvel : pushout; - // deal with bounce - if ((surface_mode & dContactBounce) != 0) { - // calculate outgoing velocity (-ve for incoming contact) - double outgoing = - dCalcVectorDot3(J1A, J1Ofs + ROW_NORMAL * rowskip + GI2__JL_MIN, node[0].body.lvel) - + dCalcVectorDot3(J1A, J1Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, node[0].body.avel); + // c1,c2 = contact points with respect to body PORs + DVector3 c1 = new DVector3(), c2 = new DVector3(); - if (b1 != null) { - outgoing += dCalcVectorDot3(J2A, J2Ofs + ROW_NORMAL * rowskip + GI2__JL_MIN, node[1].body.lvel) - + dCalcVectorDot3(J2A, J2Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, node[1].body.avel); + // get normal, with sign adjusted for body1/body2 polarity + DVector3 normal = new DVector3(); + if ((flags & dJOINT_REVERSE) != 0) { + dCopyNegatedVector3(normal, contact.geom.normal); + } + else { + dCopyVector3(normal, contact.geom.normal); + } + + DxBody b1 = node[1].body; + if (b1 != null) { + dSubtractVectors3(c2, contact.geom.pos, b1.posr().pos()); + // set Jacobian for b1 normal + dCopyNegatedVector3(J2A, J2Ofs + ROW_NORMAL * rowskip + GI2__JL_MIN, normal); + dCalcVectorCross3(J2A, J2Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, normal, c2); //== dCalcVectorCross3( J2 + GI2__JA_MIN, c2, normal ); dNegateVector3( J2 + GI2__JA_MIN ); + if (apply_bounce) { + outgoing /*+*/= dCalcVectorDot3(J2A, J2Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, node[1].body.avel) + - dCalcVectorDot3(normal, node[1].body.lvel); } + } + + DxBody b0 = node[0].body; + dSubtractVectors3(c1, contact.geom.pos, b0.posr().pos()); + // set Jacobian for b0 normal + dCopyVector3(J1A, J1Ofs + ROW_NORMAL * rowskip + GI2__JL_MIN, normal); + dCalcVectorCross3(J1A, J1Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, c1, normal); + if (apply_bounce) { + // calculate outgoing velocity (-ve for incoming contact) + outgoing += dCalcVectorDot3(J1A, J1Ofs + ROW_NORMAL * rowskip + GI2__JA_MIN, node[0].body.avel) + + dCalcVectorDot3(normal, node[0].body.lvel); + } - outgoing -= motionN; + // deal with bounce + if (apply_bounce) { + double negated_outgoing = motionN - outgoing; // only apply bounce if the outgoing velocity is greater than the // threshold, and if the resulting c[rowNormal] exceeds what we already have. - if (contact.surface.bounce_vel >= 0 && (-outgoing) > contact.surface.bounce_vel) { - final double newc = -contact.surface.bounce * outgoing + motionN; - if (newc > c) c = newc; + dIASSERT(contact.surface.bounce_vel >= 0); + if (/*contact.surface.bounce_vel >= 0 &&*/ + negated_outgoing > contact.surface.bounce_vel) { + final double newc = contact.surface.bounce * negated_outgoing + motionN; + if (newc > c) { c = newc; } } } @@ -218,7 +253,9 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A int currRowSkip = row * rowskip, currPairSkip = row * pairskip; // first friction direction - if (contact.surface.mu > 0) { + final double mu = contact.surface.mu; + + if (mu > 0) { dCopyVector3(J1A, J1Ofs + currRowSkip + GI2__JL_MIN, t1); dCalcVectorCross3(J1A, J1Ofs + currRowSkip + GI2__JA_MIN, c1, t1); @@ -239,8 +276,8 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A // set LCP bounds and friction index. this depends on the approximation // mode - pairLoHiA[pairLoHiOfs + currPairSkip + GI2_LO] = -contact.surface.mu; - pairLoHiA[pairLoHiOfs + currPairSkip + GI2_HI] = contact.surface.mu; + pairLoHiA[pairLoHiOfs + currPairSkip + GI2_LO] = -mu; + pairLoHiA[pairLoHiOfs + currPairSkip + GI2_HI] = mu; if ((surface_mode & dContactApprox1_1) != 0) { findexA[findexOfs + row] = 0; @@ -251,9 +288,8 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A currPairSkip += pairskip; } - final double mu2 = (surface_mode & dContactMu2) != 0 ? contact.surface.mu2 : contact.surface.mu; - // second friction direction + final double mu2 = (surface_mode & dContactMu2) != 0 ? contact.surface.mu2 : mu; if (mu2 > 0) { dCopyVector3(J1A, J1Ofs + currRowSkip + GI2__JL_MIN, t2); dCalcVectorCross3(J1A, J1Ofs + currRowSkip + GI2__JA_MIN, c1, t2); @@ -289,6 +325,15 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A // Handle rolling/spinning friction if ((surface_mode & dContactRolling) != 0) { + + DVector3C[] ax = { + t1, // Rolling around t1 creates movement parallel to t2 + t2, + normal // Spinning axis + }; + + final int[] approx_bits = { dContactApprox1_1, dContactApprox1_2, dContactApprox1_N }; + // Get the coefficients double[] rho = new double[3]; rho[0] = contact.surface.rho; @@ -300,17 +345,6 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A rho[2] = rho[0]; } - final DVector3C[] ax = new DVector3[3]; - ax[0] = t1; // Rolling around t1 creates movement parallel to t2 - ax[1] = t2; - ax[2] = normal; // Spinning axis - - // Should we use proportional force? - boolean[] approx = new boolean[3]; - approx[0] = (surface_mode & dContactApprox1_1) != 0; - approx[1] = (surface_mode & dContactApprox1_2) != 0; - approx[2] = (surface_mode & dContactApprox1_N) != 0; - for (int i = 0; i != 3; ++i) { if (rho[i] > 0) { // Set the angular axis @@ -324,8 +358,9 @@ public void getInfo2(double worldFPS, double worldERP, int rowskip, double[] J1A pairLoHiA[pairLoHiOfs + currPairSkip + GI2_LO] = -rho[i]; pairLoHiA[pairLoHiOfs + currPairSkip + GI2_HI] = rho[i]; - // Make limits proportional to normal force - if (approx[i]) { + // Should we use proportional force? + if ((surface_mode & approx_bits[i]) != 0) { + // Make limits proportional to normal force findexA[findexOfs + row] = 0; } diff --git a/core/src/main/java/org/ode4j/ode/internal/joints/DxJointHinge2.java b/core/src/main/java/org/ode4j/ode/internal/joints/DxJointHinge2.java index ec18a188..8687b70d 100644 --- a/core/src/main/java/org/ode4j/ode/internal/joints/DxJointHinge2.java +++ b/core/src/main/java/org/ode4j/ode/internal/joints/DxJointHinge2.java @@ -392,6 +392,12 @@ private void dJointGetHinge2Axis1( DVector3 result ) { dMultiply0_331( result, node[0].body.posr().R(), _axis1 ); } + else + { + // dZeroVector3(result); + result.setZero(); + dUASSERT( false, "the joint does not have first body attached" ); + } } @@ -402,6 +408,12 @@ private void dJointGetHinge2Axis2( DVector3 result ) { dMultiply0_331( result, node[1].body.posr().R(), _axis2 ); } + else + { + // dZeroVector3(result); + result.setZero(); + dUASSERT( false, "the joint does not have second body attached" ); + } } diff --git a/core/src/main/java/org/ode4j/ode/internal/joints/DxJointPU.java b/core/src/main/java/org/ode4j/ode/internal/joints/DxJointPU.java index 63c7bbd7..795a07e4 100644 --- a/core/src/main/java/org/ode4j/ode/internal/joints/DxJointPU.java +++ b/core/src/main/java/org/ode4j/ode/internal/joints/DxJointPU.java @@ -208,6 +208,10 @@ public double dJointGetPUPositionRate() if ( node[1].body!=null ) { // Find joint->anchor2 in global coordinates + + // NOTE! anchor2 needs a volatile assignment on the multiplication to discard computation errors. + // Otherwise, tests fail for single type on x86. + // dxTruncToType::dMultiply0_331(anchor2, joint->node[1].body->posr.R, joint->anchor2); dMultiply0_331( anchor2, node[1].body.posr().R(), _anchor2 ); // r.v[0] = ( node[0].body._posr.pos.v[0] - diff --git a/core/src/main/java/org/ode4j/ode/internal/libccd/CCD.java b/core/src/main/java/org/ode4j/ode/internal/libccd/CCD.java index feb15e65..d9415dfb 100644 --- a/core/src/main/java/org/ode4j/ode/internal/libccd/CCD.java +++ b/core/src/main/java/org/ode4j/ode/internal/libccd/CCD.java @@ -26,6 +26,7 @@ import org.ode4j.ode.internal.libccd.CCDSimplex.ccd_simplex_t; import org.ode4j.ode.internal.libccd.CCDVec3.ccd_vec3_t; +import static org.ode4j.ode.internal.libccd.CCDCustomVec3.ccdVec3SafeNormalize; import static org.ode4j.ode.internal.libccd.CCDPolyTope.*; import static org.ode4j.ode.internal.libccd.CCDSimplex.*; import static org.ode4j.ode.internal.libccd.CCDSupport.*; @@ -242,7 +243,7 @@ public static int ccdGJKPenetration(final Object obj1, final Object obj2, final RefDouble depth, ccd_vec3_t dir, ccd_vec3_t pos) { ccd_pt_t polytope = new ccd_pt_t(); - final Ref> nearestRef = new Ref>(); + final Ref> nearestRef = new Ref<>(); int ret; ccdPtInit(polytope); @@ -251,15 +252,16 @@ public static int ccdGJKPenetration(final Object obj1, final Object obj2, final // set separation vector if (ret == 0 && nearestRef.get() != null){ - // compute depth of penetration - depth.set( CCD_SQRT(nearestRef.get().dist) ); - // store normalized direction vector ccdVec3Copy(dir, nearestRef.get().witness); - ccdVec3Normalize(dir); + ret = ccdVec3SafeNormalize(dir); - // compute position - penEPAPos(polytope, nearestRef.get(), pos); + if (ret == 0) { + // compute depth of penetration + depth.d = CCD_SQRT(nearestRef.get().dist); + // compute position + penEPAPos(polytope, nearestRef.get(), pos); + } } ccdPtDestroy(polytope); diff --git a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDCustomQuat.java b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDCustomQuat.java new file mode 100644 index 00000000..50479a92 --- /dev/null +++ b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDCustomQuat.java @@ -0,0 +1,59 @@ +/*** + * libccd + * --------------------------------- + * Copyright (c)2010 Daniel Fiser + * Java-port: Copyright (c) 2009-2014 Tilmann Zaeschke + * + * + * This file is part of libccd. + * + * Distributed under the OSI-approved BSD License (the "License"); + * see accompanying file BDS-LICENSE for details or see + * . + * + * This software is distributed WITHOUT ANY WARRANTY; without even the + * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the License for more information. + */ +package org.ode4j.ode.internal.libccd; + +import static org.ode4j.ode.internal.libccd.CCDQuat.*; +import static org.ode4j.ode.internal.libccd.CCDVec3.*; + +public class CCDCustomQuat { + + /** + * Rotate vector s by quaternion q and put result into d. + */ + // _ccd_inline void ccdQuatRotVec2(ccd_vec3_t *d, const ccd_vec3_t *s, const ccd_quat_t *q); + public static void ccdQuatRotVec2(ccd_vec3_t d, final ccd_vec3_t s, final ccd_quat_t q) { +//#ifndef dLIBCCD_USE_SYSTEM +// // original version: 31 mul + 21 add +// // optimized version: 18 mul + 12 add +// // formula: d = s + 2 * cross(q.xyz, cross(q.xyz, v) + q.w * s) +// ccd_real_t cross1_x, cross1_y, cross1_z, cross2_x, cross2_y, cross2_z; +// ccd_real_t x, y, z, w; +// ccd_real_t vx, vy, vz; +// +// vx = ccdVec3X(s); +// vy = ccdVec3Y(s); +// vz = ccdVec3Z(s); +// +// w = q->q[3]; +// x = q->q[0]; +// y = q->q[1]; +// z = q->q[2]; +// +// cross1_x = y * vz - z * vy + w * vx; +// cross1_y = z * vx - x * vz + w * vy; +// cross1_z = x * vy - y * vx + w * vz; +// cross2_x = y * cross1_z - z * cross1_y; +// cross2_y = z * cross1_x - x * cross1_z; +// cross2_z = x * cross1_y - y * cross1_x; +// ccdVec3Set(d, vx + 2 * cross2_x, vy + 2 * cross2_y, vz + 2 * cross2_z); +//#else + ccdVec3Copy(d, s); + ccdQuatRotVec(d, q); +//#endif + } +} diff --git a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDCustomVec3.java b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDCustomVec3.java new file mode 100644 index 00000000..7f245fe8 --- /dev/null +++ b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDCustomVec3.java @@ -0,0 +1,86 @@ +/*** + * libccd + * --------------------------------- + * Copyright (c)2010,2011 Daniel Fiser + * Java-port: Copyright (c) 2009-2014 Tilmann Zaeschke + * + * + * This file is part of libccd. + * + * Distributed under the OSI-approved BSD License (the "License"); + * see accompanying file BDS-LICENSE for details or see + * . + * + * This software is distributed WITHOUT ANY WARRANTY; without even the + * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the License for more information. + */ +package org.ode4j.ode.internal.libccd; + +import static org.ode4j.ode.internal.libccd.CCDVec3.*; + +public class CCDCustomVec3 { + + // # define CCD_ATAN2(x, y) (atan2((x), (y))) /*!< atan2 of two floats */ + public static double CCD_ATAN2(double x, double y) { + return Math.atan2(x, y); + } + + /** + * d = v + w + */ + public static void ccdVec3Add2(ccd_vec3_t d, final ccd_vec3_t v, final ccd_vec3_t w) { +//#ifndef dLIBCCD_USE_SYSTEM +// d->v[0] = v->v[0] + w->v[0]; +// d->v[1] = v->v[1] + w->v[1]; +// d->v[2] = v->v[2] + w->v[2]; +//#else + ccdVec3Copy(d, v); + ccdVec3Add(d, w); +//#endif + } + + /** + * d = s * k; + */ + public static void ccdVec3CopyScaled(ccd_vec3_t d, final ccd_vec3_t s, double k) { +//#ifndef dLIBCCD_USE_SYSTEM +// d->v[0] = s->v[0] * k; +// d->v[1] = s->v[1] * k; +// d->v[2] = s->v[2] * k; +//#else + ccdVec3Copy(d, s); + ccdVec3Scale(d, k); +//#endif + } + + /** + * d = v + s * k; + */ + public static void ccdVec3AddScaled(ccd_vec3_t d, final ccd_vec3_t v, final ccd_vec3_t s, double k) { +//#ifndef dLIBCCD_USE_SYSTEM +// d->v[0] = v->v[0] + s->v[0] * k; +// d->v[1] = v->v[1] + s->v[1] * k; +// d->v[2] = v->v[2] + s->v[2] * k; +//#else + ccdVec3Copy(d, s); + ccdVec3Scale(d, k); + ccdVec3Add(d, v); +//#endif + } + + /** + * Normalizes given vector to unit length. + */ + public static int ccdVec3SafeNormalize(ccd_vec3_t d) { + int result = -1; + + double len = CCD_SQRT(ccdVec3Len2(d)); + if (len >= CCD_EPS) { + ccdVec3Scale(d, CCD_ONE / len); + result = 0; + } + + return result; + } +} diff --git a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDMPR.java b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDMPR.java index d6433d07..036bff51 100644 --- a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDMPR.java +++ b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDMPR.java @@ -20,6 +20,7 @@ import org.ode4j.ode.internal.cpp4j.java.RefDouble; import org.ode4j.ode.internal.libccd.CCD.ccd_t; +import static org.ode4j.ode.internal.libccd.CCDCustomVec3.ccdVec3SafeNormalize; import static org.ode4j.ode.internal.libccd.CCDSimplex.*; import static org.ode4j.ode.internal.libccd.CCDSupport.*; import static org.ode4j.ode.internal.libccd.CCDVec3.*; @@ -86,22 +87,27 @@ public static int ccdMPRPenetration(final Object obj1, final Object obj2, final // Origin isn't inside portal - no collision. return -1; - }else if (res == 1){ + } else if (res == 1) { // Touching contact on portal's v1. findPenetrTouch(obj1, obj2, ccd, portal, depth, dir, pos); - }else if (res == 2){ + } else if (res == 2) { // Origin lies on v0-v1 segment. - findPenetrSegment(obj1, obj2, ccd, portal, depth, dir, pos); + if (findPenetrSegment(obj1, obj2, ccd, portal, depth, dir, pos) != 0) { + return -1; + } - }else if (res == 0){ + } else if (res == 0) { // Phase 2: Portal refinement res = refinePortal(obj1, obj2, ccd, portal); - if (res < 0) - return -1; + if (res < 0) { + return -1; + } // Phase 3. Penetration info - findPenetr(obj1, obj2, ccd, portal, depth, dir, pos); + if (findPenetr(obj1, obj2, ccd, portal, depth, dir, pos) != 0) { + return -1; + } } return 0; @@ -152,8 +158,10 @@ private static int discoverPortal(final Object obj1, final Object obj2, // vertex 1 = support in direction of origin ccdVec3Copy(dir, ccdSimplexPoint0(portal).v); ccdVec3Scale(dir, (-1.)); - ccdVec3Normalize(dir); - __ccdSupport(obj1, obj2, dir, ccd, ccdSimplexPointW1(portal)); + if (ccdVec3SafeNormalize(dir) != 0) { + return -1; + } + __ccdSupport(obj1, obj2, dir, ccd, ccdSimplexPointW1(portal)); ccdSimplexSetSize(portal,2); // test if origin isn't outside of v1 @@ -175,11 +183,14 @@ private static int discoverPortal(final Object obj1, final Object obj2, } } - ccdVec3Normalize(dir); + if (ccdVec3SafeNormalize(dir) != 0) { + return -1; + } __ccdSupport(obj1, obj2, dir, ccd, ccdSimplexPointW2(portal)); dot = ccdVec3Dot(ccdSimplexPoint2(portal).v, dir); - if (ccdIsZero(dot) || dot < CCD_ZERO) - return -1; + if (ccdIsZero(dot) || dot < CCD_ZERO) { + return -1; + } ccdSimplexSetSize(portal, 3); @@ -189,7 +200,9 @@ private static int discoverPortal(final Object obj1, final Object obj2, ccdVec3Sub2(vb, ccdSimplexPoint2(portal).v, ccdSimplexPoint0(portal).v); ccdVec3Cross(dir, va, vb); - ccdVec3Normalize(dir); + if (ccdVec3SafeNormalize(dir) != 0) { + return -1; + } // it is better to form portal faces to be oriented "outside" origin dot = ccdVec3Dot(dir, ccdSimplexPoint0(portal).v); @@ -201,8 +214,9 @@ private static int discoverPortal(final Object obj1, final Object obj2, while (ccdSimplexSize(portal) < 4){ __ccdSupport(obj1, obj2, dir, ccd, ccdSimplexPointW3(portal)); dot = ccdVec3Dot(ccdSimplexPoint3(portal).v, dir); - if (ccdIsZero(dot) || dot < CCD_ZERO) - return -1; + if (ccdIsZero(dot) || dot < CCD_ZERO) { + return -1; + } cont = 0; @@ -234,8 +248,10 @@ private static int discoverPortal(final Object obj1, final Object obj2, ccdVec3Sub2(vb, ccdSimplexPoint2(portal).v, ccdSimplexPoint0(portal).v); ccdVec3Cross(dir, va, vb); - ccdVec3Normalize(dir); - }else{ + if (ccdVec3SafeNormalize(dir) != 0) { + return -1; + } + } else { ccdSimplexSetSize(portal, 4); } } @@ -255,7 +271,9 @@ private static int refinePortal(final Object obj1, final Object obj2, while (true){ // compute direction outside the portal (from v0 throught v1,v2,v3 // face) - portalDir(portal, dir); + if (portalDir(portal, dir) != 0) { + return -1; + } // test if origin is inside the portal if (portalEncapsulesOrigin(portal, dir)) @@ -282,7 +300,7 @@ private static int refinePortal(final Object obj1, final Object obj2, /** Finds penetration info by expanding provided portal. */ - private static void findPenetr(final Object obj1, final Object obj2, final ccd_t ccd, + private static int findPenetr(final Object obj1, final Object obj2, final ccd_t ccd, ccd_simplex_t portal, RefDouble depth, ccd_vec3_t pdir, ccd_vec3_t pos) { @@ -293,8 +311,11 @@ private static void findPenetr(final Object obj1, final Object obj2, final ccd_t iterations = 0L; while (true){ // compute portal direction and obtain next support point - portalDir(portal, dir); - __ccdSupport(obj1, obj2, dir, ccd, v4); + if (portalDir(portal, dir) != 0) { + return -1; + } + + __ccdSupport(obj1, obj2, dir, ccd, v4); // reached tolerance . find penetration info if (portalReachTolerance(portal, v4, dir, ccd) @@ -305,12 +326,16 @@ private static void findPenetr(final Object obj1, final Object obj2, final ccd_t ccdSimplexPoint3(portal).v, pdir)); depth.set( CCD_SQRT(depth.get())); - ccdVec3Normalize(pdir); + if (ccdVec3SafeNormalize(pdir) != 0) { + return -1; + } // barycentric coordinates: - findPos(obj1, obj2, ccd, portal, pos); + if (findPos(obj1, obj2, ccd, portal, pos) != 0) { + return -1; + } - return; + return 0; } expandPortal(portal, v4); @@ -336,7 +361,7 @@ private static void findPenetrTouch(final Object obj1, final Object obj2, final } /** Find penetration info if origin lies on portal's segment v0-v1 */ - private static void findPenetrSegment(final Object obj1, final Object obj2, final ccd_t ccd, + private static int findPenetrSegment(final Object obj1, final Object obj2, final ccd_t ccd, ccd_simplex_t portal, RefDouble depth, ccd_vec3_t dir, ccd_vec3_t pos) { @@ -364,13 +389,16 @@ private static void findPenetrSegment(final Object obj1, final Object obj2, fina ccdVec3Copy(dir, ccdSimplexPoint1(portal).v); depth.set( CCD_SQRT(ccdVec3Len2(dir)) ); - ccdVec3Normalize(dir); + if (ccdVec3SafeNormalize(dir) != 0) { + return -1; + } + return 0; } /** Finds position vector from fully established portal */ - private static void findPos(final Object obj1, final Object obj2, final ccd_t ccd, + private static int findPos(final Object obj1, final Object obj2, final ccd_t ccd, final ccd_simplex_t portal, ccd_vec3_t pos) { final ccd_vec3_t dir = new ccd_vec3_t(); @@ -379,7 +407,9 @@ private static void findPos(final Object obj1, final Object obj2, final ccd_t cc double sum, inv; final ccd_vec3_t vec = new ccd_vec3_t(), p1 = new ccd_vec3_t(), p2 = new ccd_vec3_t(); - portalDir(portal, dir); + if (portalDir(portal, dir) != 0) { + return -1; + } // use barycentric coordinates of tetrahedron to find origin ccdVec3Cross(vec, ccdSimplexPoint1(portal).v, @@ -458,12 +488,13 @@ private static void findPos(final Object obj1, final Object obj2, final ccd_t cc ccdVec3Copy(pos, p1); ccdVec3Add(pos, p2); ccdVec3Scale(pos, 0.5); + return 0; } /** Extends portal with new support point. * Portal must have face v1-v2-v3 arranged to face outside portal. */ - private static final void expandPortal(ccd_simplex_t portal, + private static void expandPortal(ccd_simplex_t portal, final ccd_support_t v4) { double dot; @@ -491,7 +522,7 @@ private static final void expandPortal(ccd_simplex_t portal, /** Fill dir with direction outside portal. Portal's v1-v2-v3 face must be * arranged in correct order! */ - private static final void portalDir(final ccd_simplex_t portal, ccd_vec3_t dir) + private static int portalDir(final ccd_simplex_t portal, ccd_vec3_t dir) { final ccd_vec3_t v2v1 = new ccd_vec3_t(), v3v1 = new ccd_vec3_t(); @@ -500,7 +531,10 @@ private static final void portalDir(final ccd_simplex_t portal, ccd_vec3_t dir) ccdVec3Sub2(v3v1, ccdSimplexPoint3(portal).v, ccdSimplexPoint1(portal).v); ccdVec3Cross(dir, v2v1, v3v1); - ccdVec3Normalize(dir); + if (ccdVec3SafeNormalize(dir) != 0) { + return -1; + } + return 0; } /** Returns true if portal encapsules origin (0,0,0), dir is direction of diff --git a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDQuat.java b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDQuat.java index 00e415a6..b02f3c4a 100644 --- a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDQuat.java +++ b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDQuat.java @@ -198,19 +198,19 @@ public static final int ccdQuatInvert2(ccd_quat_t dest, final ccd_quat_t src) } /** - * Rotate vector s by quaternion q and put result into d. + * Rotate vector v by quaternion q. */ - public static void ccdQuatRotVec2(ccd_vec3_t d, final ccd_vec3_t s, final ccd_quat_t q) { + public static void ccdQuatRotVec(ccd_vec3_t v, final ccd_quat_t q) { // original version: 31 mul + 21 add // optimized version: 18 mul + 12 add - // formula: d = s + 2 * cross(q.xyz, cross(q.xyz, v) + q.w * s) + // formula: v = v + 2 * cross(q.xyz, cross(q.xyz, v) + q.w * v) double cross1_x, cross1_y, cross1_z, cross2_x, cross2_y, cross2_z; double x, y, z, w; double vx, vy, vz; - vx = ccdVec3X(s); - vy = ccdVec3Y(s); - vz = ccdVec3Z(s); + vx = ccdVec3X(v); + vy = ccdVec3Y(v); + vz = ccdVec3Z(v); w = q.q3; x = q.q0; @@ -223,14 +223,7 @@ public static void ccdQuatRotVec2(ccd_vec3_t d, final ccd_vec3_t s, final ccd_qu cross2_x = y * cross1_z - z * cross1_y; cross2_y = z * cross1_x - x * cross1_z; cross2_z = x * cross1_y - y * cross1_x; - ccdVec3Set(d, vx + 2 * cross2_x, vy + 2 * cross2_y, vz + 2 * cross2_z); - } - - /** - * Rotate vector v by quaternion q. - */ - public static void ccdQuatRotVec(ccd_vec3_t v, final ccd_quat_t q) { - ccdQuatRotVec2(v, v, q); + ccdVec3Set(v, vx + 2 * cross2_x, vy + 2 * cross2_y, vz + 2 * cross2_z); } } diff --git a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDVec3.java b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDVec3.java index 22f98452..2cbbe656 100644 --- a/core/src/main/java/org/ode4j/ode/internal/libccd/CCDVec3.java +++ b/core/src/main/java/org/ode4j/ode/internal/libccd/CCDVec3.java @@ -39,10 +39,6 @@ public class CCDVec3 { /* minimum of two floats */ static double CCD_FMIN(double x, double y) { return Math.min(x, y); } - public static double CCD_ATAN2(double x, double y) { - return Math.atan2(x, y); - } - //#define CCD_ONE CCD_REAL(1.) public static final double CCD_ONE = 1; //#define CCD_ZERO CCD_REAL(0.) @@ -262,15 +258,6 @@ public static void ccdVec3Sub2(ccd_vec3_t d, final ccd_vec3_t v, final ccd_vec3_ d.v2 = v.v2 - w.v2; } - /** - * d = v + w - */ - public static void ccdVec3Add2(ccd_vec3_t d, final ccd_vec3_t v, final ccd_vec3_t w) { - d.v0 = v.v0 + w.v0; - d.v1 = v.v1 + w.v1; - d.v2 = v.v2 + w.v2; - } - /** * d = d * k; * @@ -283,40 +270,13 @@ public static void ccdVec3Scale(ccd_vec3_t d, double k) { d.v2 *= k; } - /** - * d = s * k; - */ - public static void ccdVec3CopyScaled(ccd_vec3_t d, final ccd_vec3_t s, double k) { - d.v0 = s.v0 * k; - d.v1 = s.v1 * k; - d.v2 = s.v2 * k; - } - - /** - * d = v + s * k; - */ - static void ccdVec3AddScaled(ccd_vec3_t d, final ccd_vec3_t v, final ccd_vec3_t s, double k) { - d.v0 = v.v0 + s.v0 * k; - d.v1 = v.v1 + s.v1 * k; - d.v2 = v.v2 + s.v2 * k; - } - /** * Normalizes given vector to unit length. - * - * WARNING: This returns TRUE for a failure! - * - * @param d d - * @return 'true' if normalization failed. */ - public static boolean ccdVec3Normalize(ccd_vec3_t d) + public static void ccdVec3Normalize(ccd_vec3_t d) { - double len = CCD_SQRT(ccdVec3Len2(d)); - if (len > CCD_EPS) { - ccdVec3Scale(d, CCD_ONE / len); - return false; - } - return true; + double k = CCD_ONE / CCD_SQRT(ccdVec3Len2(d)); + ccdVec3Scale(d, k); } /** diff --git a/core/src/main/java/org/ode4j/ode/internal/processmem/DxUtil.java b/core/src/main/java/org/ode4j/ode/internal/processmem/DxUtil.java index 2d050f3a..1fc067b2 100644 --- a/core/src/main/java/org/ode4j/ode/internal/processmem/DxUtil.java +++ b/core/src/main/java/org/ode4j/ode/internal/processmem/DxUtil.java @@ -58,23 +58,22 @@ public class DxUtil { 1.2f, //dWORLDSTEP_RESERVEFACTOR_DEFAULT, 65536);// dWORLDSTEP_RESERVESIZE_DEFAULT); - - static final int sizeof(Class cls) { - if (cls == DxWorldProcessMemArena.class) { - return -1; - } - return -1; - } - - public static final int sizeof(Object o) { - if (o instanceof Class) { - return sizeof((Class)o); - } - return sizeof(o.getClass()); - } - + + /* the efficient alignment. most platforms align data structures to some + * number of bytes, but this is not always the most efficient alignment. + * for example, many x86 compilers align to 4 bytes, but on a pentium it + * is important to align doubles to 8 byte boundaries (for speed), and + * the 4 floats in a SIMD register to 16 byte boundaries. many other + * platforms have similar behavior. setting a larger alignment can waste + * a (very) small amount of memory. NOTE: this number must be a power of + * two. this is set to 16 by default. + */ public static final int EFFICIENT_ALIGNMENT = 16; - + + /* ********************************************************************************************************* */ + /* ************************* TZ: This has been move here from Common because it is mostly required here **** */ + /* ********************************************************************************************************* */ + // #define dEFFICIENT_SIZE(x) (((x)+(EFFICIENT_ALIGNMENT-1)) & ~((size_t)(EFFICIENT_ALIGNMENT-1))) static final int dEFFICIENT_SIZE(int x) { return x; diff --git a/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessContext.java b/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessContext.java index 794a637d..39653abe 100644 --- a/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessContext.java +++ b/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessContext.java @@ -105,6 +105,14 @@ public void DESTRUCTOR() { Common.dIASSERT((m_pswObjectsAllocWorld != null) == (m_pmgStepperMutexGroup != null)); + // TODO TZ: Ignored for migration to 0.16.2 + // if (m_pswObjectsAllocWorld != null) + // { + // throw new UnsupportedOperationException("This has not been ported yet from ODE"); + // // TODO m_pswObjectsAllocWorld->FreeMutexGroup(m_pmgStepperMutexGroup); + // // m_pswObjectsAllocWorld->FreeThreadedCallWait(m_pcwIslandsSteppingWait); -- The stock call wait can not be freed + // } + DxWorldProcessMemArena pmaStepperArenas = m_pmaStepperArenas.get(); if (pmaStepperArenas != null) { @@ -396,14 +404,16 @@ public static boolean dxReallocateWorldProcessContext (DxWorld world, DxWorldPro //TODO this is just wrong!!!! (TZ) int stepperReqWithCallContext = stepperReq + DxUtil.dEFFICIENT_SIZE(1);//DxSingleIslandCallContext.class); - int islandThreadsCount = world.GetThreadingIslandsMaxThreadsCount(null); - if (!context.ReallocateStepperMemArenas(world, islandThreadsCount, stepperReqWithCallContext, - memmgr, reserveInfo.m_fReserveFactor, reserveInfo.m_uiReserveMinimum)) - { - break; - } + // TODO TZ: This has not been migrfated from ODE 0.16.2 + // throw new UnsupportedOperationException("This has not been ported yet from ODE"); + int islandThreadsCount = world.GetThreadingIslandsMaxThreadsCount(null); + if (!context.ReallocateStepperMemArenas(world, islandThreadsCount, stepperReqWithCallContext, + memmgr, reserveInfo.m_fReserveFactor, reserveInfo.m_uiReserveMinimum)) + { + break; + } - result = true; + result = true; } while (false); diff --git a/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessMemArena.java b/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessMemArena.java index 7ab6ce57..2de04472 100644 --- a/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessMemArena.java +++ b/core/src/main/java/org/ode4j/ode/internal/processmem/DxWorldProcessMemArena.java @@ -36,8 +36,8 @@ public final class DxWorldProcessMemArena { // public: // #define BUFFER_TO_ARENA_EXTRA (EFFICIENT_ALIGNMENT + dEFFICIENT_SIZE(sizeof(dxWorldProcessMemArena))) private final static int BUFFER_TO_ARENA_EXTRA () { - return (DxUtil.EFFICIENT_ALIGNMENT + DxUtil.dEFFICIENT_SIZE( - DxUtil.sizeof(DxWorldProcessMemArena.class))); + return 0; //(DxUtil.EFFICIENT_ALIGNMENT + DxUtil.dEFFICIENT_SIZE( + //DxUtil.sizeof(DxWorldProcessMemArena.class))); } static boolean IsArenaPossible(int nBufferSize) diff --git a/core/src/main/java/org/ode4j/ode/threading/Atomics.java b/core/src/main/java/org/ode4j/ode/threading/Atomics.java index 759ee060..7206fb02 100644 --- a/core/src/main/java/org/ode4j/ode/threading/Atomics.java +++ b/core/src/main/java/org/ode4j/ode/threading/Atomics.java @@ -57,7 +57,10 @@ public static int ThrsafeIncrementIntUpToLimit(AtomicInteger storagePointer, int int resultValue; while (true) { resultValue = storagePointer.get(); - if (resultValue == limitValue) { + // The ">=" comparison is used here to allow continuing incrementing the destination + // without waiting for all the threads to pass the barrier of checking its value + if (resultValue >= limitValue) { + resultValue = limitValue; break; } if (ThrsafeCompareExchange(storagePointer, resultValue, resultValue + 1)) { @@ -72,7 +75,10 @@ public static int ThrsafeIncrementSizeUpToLimit(AtomicInteger storagePointer, in int resultValue; while (true) { resultValue = storagePointer.get(); - if (resultValue == limitValue) { + // The ">=" comparison is not required here at present ("==" could be used). + // It is just used this way to match the other function above. + if (resultValue >= limitValue) { + resultValue = limitValue; break; } //if (ThrsafeCompareExchangePointer((volatile atomicptr *)storagePointer, (atomicptr)resultValue, (atomicptr)(resultValue + 1))) { diff --git a/demo-cpp/src/test/java/org/ode4j/tests/CollisionTest.java b/demo-cpp/src/test/java/org/ode4j/tests/CollisionTest.java index f3b68460..fcd38f62 100644 --- a/demo-cpp/src/test/java/org/ode4j/tests/CollisionTest.java +++ b/demo-cpp/src/test/java/org/ode4j/tests/CollisionTest.java @@ -24,33 +24,18 @@ *************************************************************************/ package org.ode4j.tests; -import static org.ode4j.cpp.internal.ApiCppCollision.dCollide; -import static org.ode4j.cpp.internal.ApiCppCollision.dCreateHeightfield; -import static org.ode4j.cpp.internal.ApiCppCollision.dCreateRay; -import static org.ode4j.cpp.internal.ApiCppCollision.dCreateSphere; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomDestroy; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomHeightfieldDataBuildByte; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomHeightfieldDataCreate; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomHeightfieldDataDestroy; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomHeightfieldDataSetBounds; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomRaySet; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomSetPosition; -import static org.ode4j.cpp.internal.ApiCppCollision.dGeomSetRotation; +import static org.ode4j.cpp.internal.ApiCppCollision.*; import static org.ode4j.cpp.internal.ApiCppCollisionTrimesh.dCreateTriMesh; import static org.ode4j.cpp.internal.ApiCppCollisionTrimesh.dGeomTriMeshDataBuildSingle; import static org.ode4j.cpp.internal.ApiCppCollisionTrimesh.dGeomTriMeshDataCreate; -import static org.ode4j.tests.UnitTestPlusPlus.CheckMacros.CHECK_ARRAY_EQUAL; -import static org.ode4j.tests.UnitTestPlusPlus.CheckMacros.CHECK_EQUAL; +import static org.ode4j.demo.ConvexPrism.*; +import static org.ode4j.ode.internal.Common.*; +import static org.ode4j.tests.UnitTestPlusPlus.CheckMacros.*; import org.junit.Test; import org.ode4j.math.DMatrix3; import org.ode4j.math.DVector3; -import org.ode4j.ode.DContactBuffer; -import org.ode4j.ode.DContactGeomBuffer; -import org.ode4j.ode.DGeom; -import org.ode4j.ode.DHeightfieldData; -import org.ode4j.ode.DRay; -import org.ode4j.ode.DTriMeshData; +import org.ode4j.ode.*; import org.ode4j.tests.UnitTestPlusPlus.TestSuperClass; /** @@ -198,5 +183,102 @@ public class CollisionTest extends TestSuperClass { dGeomHeightfieldDataDestroy(heightfieldData); } + @Test + public void test_collision_ray_convex() + { + /* + * Issue 55: ray vs convex collider does not consider the position of the convex geometry. + */ + { + DContactBuffer all_contacts = new DContactBuffer(1); + DContactGeomBuffer contacts = all_contacts.getGeomBuffer(); + DContact contact = all_contacts.get(0); // TODO correct? + + // Create convex + DGeom convex = dCreateConvex(null, + prism_planes, + prism_planecount, + prism_points, + prism_pointcount, + prism_polygons); + dGeomSetPosition(convex,0,0,0); + + // Create ray + DRay ray = dCreateRay(null, 20); + + dGeomRaySet(ray, 0, -10, 0, 0, 1, 0); + + int count = dCollide(ray, convex, 1, contacts); + + CHECK_EQUAL(1,count); + CHECK_CLOSE(0.0,contact.geom.pos.get0(), dEpsilon); + CHECK_CLOSE(-1.0,contact.geom.pos.get1(), dEpsilon); + CHECK_CLOSE(0.0,contact.geom.pos.get2(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get0(), dEpsilon); + CHECK_CLOSE(-1.0, contact.geom.normal.get1(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get2(), dEpsilon); + CHECK_CLOSE(9.0, contact.geom.depth, dEpsilon); + + // Move Ray + dGeomRaySet(ray, 5, -10, 0, 0, 1, 0); + + count = dCollide(ray, convex, 1, contacts); + + CHECK_EQUAL(1,count); + CHECK_CLOSE(5.0, contact.geom.pos.get0(), dEpsilon); + CHECK_CLOSE(-1.0, contact.geom.pos.get1(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.pos.get2(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get0(), dEpsilon); + CHECK_CLOSE(-1.0, contact.geom.normal.get1(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get2(), dEpsilon); + CHECK_CLOSE(9.0, contact.geom.depth, dEpsilon); + + // Rotate Convex + DMatrix3 rotate90z = new DMatrix3( + + 0,-1,0,0, + 1,0,0,0, + 0,0,1,0 + ); + dGeomSetRotation(convex, rotate90z); + + count = dCollide(ray, convex, 1, contacts); + + CHECK_EQUAL(0,count); + + // Move Ray + dGeomRaySet(ray, 10, 0, 0, -1, 0, 0); + count = dCollide(ray, convex, 1, contacts); + + CHECK_EQUAL(1,count); + CHECK_CLOSE(1.0, contact.geom.pos.get0(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.pos.get1(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.pos.get2(), dEpsilon); + CHECK_CLOSE(1.0, contact.geom.normal.get0(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get1(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get2(), dEpsilon); + CHECK_CLOSE(9.0,contact.geom.depth, dEpsilon); + + + // Move Ray + dGeomRaySet(ray, 10, 1000, 1000, -1, 0, 0); + // Move Geom + dGeomSetPosition(convex, 0, 1000, 1000); + + count = dCollide(ray, convex, 1, contacts); + + CHECK_EQUAL(1, count); + CHECK_CLOSE(1.0, contact.geom.pos.get0(), dEpsilon); + CHECK_CLOSE(1000.0, contact.geom.pos.get1(), dEpsilon); + CHECK_CLOSE(1000.0, contact.geom.pos.get2(), dEpsilon); + CHECK_CLOSE(1.0, contact.geom.normal.get0(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get1(), dEpsilon); + CHECK_CLOSE(0.0, contact.geom.normal.get2(), dEpsilon); + CHECK_CLOSE(9.0, contact.geom.depth, dEpsilon); + + dGeomDestroy(convex); + dGeomDestroy(ray); + } + } } diff --git a/demo-cpp/src/test/java/org/ode4j/tests/TestOdeMath.java b/demo-cpp/src/test/java/org/ode4j/tests/TestOdeMath.java index a0d13139..c16fd64e 100644 --- a/demo-cpp/src/test/java/org/ode4j/tests/TestOdeMath.java +++ b/demo-cpp/src/test/java/org/ode4j/tests/TestOdeMath.java @@ -145,8 +145,7 @@ public class TestOdeMath { v3.set( 9999999999.0, 9999.0, 9.0 ); v3.normalize(); - CHECK_EQUAL(v3.length(), 1.0); - + CHECK_CLOSE(v3.length(), 1.0, 0.001); } diff --git a/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench.java b/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench.java index 50cb3623..48a09129 100644 --- a/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench.java +++ b/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench.java @@ -273,7 +273,7 @@ public static void main(String[] args) { cycles = Integer.parseInt(args[0]); } -// fprintf(stdout, "Cycles: %u\n", cycles); +// fprintf(stdout, "Cycles: %zu\n", cycles); // fprintf(stdout, "\n"); boxbox(); diff --git a/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench2.java b/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench2.java index c980fe81..3e5d28cd 100644 --- a/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench2.java +++ b/demo-cpp/src/test/java/org/ode4j/tests/libccd/CCDTestBench2.java @@ -285,7 +285,7 @@ public static void main(String[] args) { cycles = Integer.parseInt(args[0]); } -// fprintf(stdout, "Cycles: %u\n", cycles); +// fprintf(stdout, "Cycles: %zu\n", cycles); // fprintf(stdout, "\n"); boxbox(); diff --git a/demo-cpp/src/test/java/org/ode4j/tests/math/OdeMathTZ.java b/demo-cpp/src/test/java/org/ode4j/tests/math/OdeMathTZ.java index d79689d3..4e6b0b92 100644 --- a/demo-cpp/src/test/java/org/ode4j/tests/math/OdeMathTZ.java +++ b/demo-cpp/src/test/java/org/ode4j/tests/math/OdeMathTZ.java @@ -154,8 +154,7 @@ public class OdeMathTZ { v3.set(9999999999.0, 9999.0, 9.0); v3.normalize(); - CHECK_EQUAL(v3.length(), 1.0); - + CHECK_CLOSE(v3.length(), 1.0, 0.001); } @Test public void test_dOrthogonalizeR() { diff --git a/demo/src/main/java/org/ode4j/demo/ConvexCubeGeom.java b/demo/src/main/java/org/ode4j/demo/ConvexCubeGeom.java index deaa5ecd..838a98c4 100644 --- a/demo/src/main/java/org/ode4j/demo/ConvexCubeGeom.java +++ b/demo/src/main/java/org/ode4j/demo/ConvexCubeGeom.java @@ -1,3 +1,27 @@ +/************************************************************************* + * * + * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith. * + * All rights reserved. Email: russ@q12.org Web: www.q12.org * + * Open Dynamics Engine 4J, Copyright (C) 2009-2014 Tilmann Zaeschke * + * All rights reserved. Email: ode4j@gmx.de Web: www.ode4j.org * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of EITHER: * + * (1) The GNU Lesser General Public License as published by the Free * + * Software Foundation; either version 2.1 of the License, or (at * + * your option) any later version. The text of the GNU Lesser * + * General Public License is included with this library in the * + * file LICENSE.TXT. * + * (2) The BSD-style license that is included with this library in * + * the file ODE-LICENSE-BSD.TXT and ODE4J-LICENSE-BSD.TXT. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files * + * LICENSE.TXT, ODE-LICENSE-BSD.TXT and ODE4J-LICENSE-BSD.TXT for more * + * details. * + * * + *************************************************************************/ package org.ode4j.demo; class ConvexCubeGeom { diff --git a/demo/src/main/java/org/ode4j/demo/ConvexPrism.java b/demo/src/main/java/org/ode4j/demo/ConvexPrism.java new file mode 100644 index 00000000..181f0d62 --- /dev/null +++ b/demo/src/main/java/org/ode4j/demo/ConvexPrism.java @@ -0,0 +1,60 @@ +/************************************************************************* + * * + * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith. * + * All rights reserved. Email: russ@q12.org Web: www.q12.org * + * Open Dynamics Engine 4J, Copyright (C) 2009-2014 Tilmann Zaeschke * + * All rights reserved. Email: ode4j@gmx.de Web: www.ode4j.org * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of EITHER: * + * (1) The GNU Lesser General Public License as published by the Free * + * Software Foundation; either version 2.1 of the License, or (at * + * your option) any later version. The text of the GNU Lesser * + * General Public License is included with this library in the * + * file LICENSE.TXT. * + * (2) The BSD-style license that is included with this library in * + * the file ODE-LICENSE-BSD.TXT and ODE4J-LICENSE-BSD.TXT. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files * + * LICENSE.TXT, ODE-LICENSE-BSD.TXT and ODE4J-LICENSE-BSD.TXT for more * + * details. * + * * + *************************************************************************/ +package org.ode4j.demo; + + +public class ConvexPrism { + + public static int prism_pointcount = 8; + + public static int prism_planecount = 6; + public static double[] prism_points = { // [24]= { + 10.0, 1.0, -1.0, + 10.0, -1.0, -1.0, + -10.0, -1.0, -1.0, + -10.0, 1.0, -1.0, + 10.0, 1.0, 1.0, + 10.0, -1.0, 1.0, + -10.0, -1.0, 1.0, + -10.0, 1.0, 1.0 + }; + + public static int[] prism_polygons = { + 4, 0, 1, 2, 3, + 4, 4, 7, 6, 5, + 4, 0, 4, 5, 1, + 4, 1, 5, 6, 2, + 4, 2, 6, 7, 3, + 4, 4, 0, 3, 7, + }; + public static double[] prism_planes = { + 0.0, 0.0, -1.0, 1.0, + 0.0, 0.0, 1.0, 1.0, + 1.0, 0.0, 0.0, 10.0, + 0.0, -1.0, 0.0, 1.0, + -1.0, 0.0, -0.0, 10.0, + 0.0, 1.0, 0.0, 1.0, + }; +} \ No newline at end of file diff --git a/demo/src/main/java/org/ode4j/demo/DemoCollision.java b/demo/src/main/java/org/ode4j/demo/DemoCollision.java index 4f137bbe..4fb97615 100644 --- a/demo/src/main/java/org/ode4j/demo/DemoCollision.java +++ b/demo/src/main/java/org/ode4j/demo/DemoCollision.java @@ -1105,34 +1105,59 @@ private static boolean edgeIntersectsRect (DVector3 v1, DVector3 v2, DVector3 p1, DVector3 p2, DVector3 p3) { DVector3 u1=new DVector3(),u2=new DVector3(),n=new DVector3(),tmp=new DVector3(); + //for (k=0; k<3; k++) u1[k] = p3[k]-p1[k]; u1.eqDiff(p3, p1); //for (k=0; k<3; k++) u2[k] = p2[k]-p1[k]; u2.eqDiff(p2, p1); + double d1 = dSqrt(u1.dot(u1)); double d2 = dSqrt(u2.dot(u2)); u1.normalize(); u2.normalize(); - if (dFabs(u1.dot(u2)) > 1e-6) dDebug (0,"bad u1/u2"); - n.eqCross(u1, u2); - //for (k=0; k<3; k++) tmp[k] = v2[k]-v1[k]; - tmp.eqDiff(v2, v1); - double d = -n.dot(p1); - if (dFabs(n.dot(p1)+d) > 1e-8) dDebug (0,"bad n wrt p1"); - if (dFabs(n.dot(p2)+d) > 1e-8) dDebug (0,"bad n wrt p2"); - if (dFabs(n.dot(p3)+d) > 1e-8) dDebug (0,"bad n wrt p3"); - double alpha = -(d+n.dot(v1))/n.dot(tmp); - //for (k=0; k<3; k++) tmp[k] = v1[k]+alpha*(v2[k]-v1[k]); + + double error; + // #ifdef dSINGLE + // const dReal uEpsilon = 1e-5, pEpsilon = 1e-6, tmpEpsilon = 1.5e-4; + // #else + final double uEpsilon = 1e-6, pEpsilon = 1e-8, tmpEpsilon = 1e-6; + // #endif + + error = dFabs(dCalcVectorDot3(u1, u2)); + if (error > uEpsilon) dDebug(0, "bad u1/u2"); + + dCalcVectorCross3(n, u1, u2); + + // for (k=0; k < 3; k++) tmp[k] = v2[k] - v1[k]; tmp.eqDiff(v2, v1); - tmp.eqSum(v1, tmp.scale(alpha)); - if (dFabs(n.dot(tmp)+d) > 1e-6) dDebug (0,"bad tmp"); + + double d = -dCalcVectorDot3(n, p1); + + error = dFabs(dCalcVectorDot3(n, p1) + d); + if (error > pEpsilon) dDebug(0, "bad n wrt p1"); + + error = dFabs(dCalcVectorDot3(n, p2) + d); + if (error > pEpsilon) dDebug(0, "bad n wrt p2"); + + error = dFabs(dCalcVectorDot3(n, p3) + d); + if (error > pEpsilon) dDebug(0, "bad n wrt p3"); + + double alpha = -(d + dCalcVectorDot3(n, v1)) / dCalcVectorDot3(n, tmp); + // for (k=0; k < 3; k++) tmp[k] = v1[k] + alpha * (v2[k] - v1[k]); + tmp.eqDiff(v2, v1).scale(alpha).add(v1); + + error = dFabs(dCalcVectorDot3(n, tmp) + d); + if (error > tmpEpsilon) dDebug(0, "bad tmp"); + if (alpha < 0) return false; if (alpha > 1) return false; + //for (k=0; k<3; k++) tmp[k] -= p1[k]; tmp.sub(p1); double a1 = u1.dot(tmp); double a2 = u2.dot(tmp); if (a1<0 || a2<0 || a1>d1 || a2>d2) return false; + return true; } diff --git a/demo/src/main/java/org/ode4j/demo/DemoConvex.java b/demo/src/main/java/org/ode4j/demo/DemoConvex.java new file mode 100644 index 00000000..1a2a5dbb --- /dev/null +++ b/demo/src/main/java/org/ode4j/demo/DemoConvex.java @@ -0,0 +1,312 @@ +/************************************************************************* + * * + * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith. * + * All rights reserved. Email: russ@q12.org Web: www.q12.org * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of EITHER: * + * (1) The GNU Lesser General Public License as published by the Free * + * Software Foundation; either version 2.1 of the License, or (at * + * your option) any later version. The text of the GNU Lesser * + * General Public License is included with this library in the * + * file LICENSE.TXT. * + * (2) The BSD-style license that is included with this library in * + * the file LICENSE-BSD.TXT. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files * + * LICENSE.TXT and LICENSE-BSD.TXT for more details. * + * * + *************************************************************************/ +package org.ode4j.demo; + +import org.ode4j.drawstuff.DrawStuff; +import org.ode4j.math.DMatrix3C; +import org.ode4j.math.DQuaternion; +import org.ode4j.math.DVector3; +import org.ode4j.math.DVector3C; +import org.ode4j.ode.*; + +import static org.ode4j.drawstuff.DrawStuff.*; +import static org.ode4j.drawstuff.DrawStuff.DS_TEXTURE_NUMBER.DS_NONE; +import static org.ode4j.ode.OdeConstants.*; +import static org.ode4j.ode.internal.Common.M_PI; +import static org.ode4j.demo.Halton235Geom.*; + +/** + * Convex demo. + * Serves as a test for the convex geometry. + * By Bram Stolk. + */ +public class DemoConvex extends DrawStuff.dsFunctions { + + // Height at which we drop the composite block. + private static final double H = 4.20; + + private static DWorld world; + private static DSpace space; + + private static DBody mbody; + + private static final DBody[] hbody = new DBody[halton_numc ]; + private static final DGeom[] hgeom = new DGeom[halton_numc ]; + + private static DJointGroup contactgroup; + + private static boolean drawpos = false; + private static boolean solidkernel = false; + + + private final DGeom.DNearCallback nearCallback = new DGeom.DNearCallback() { + @Override + public void call(Object data, DGeom o1, DGeom o2) { + nearCallback(data, o1, o2); + } + }; + + // this is called by dSpaceCollide when two objects in space are + // potentially colliding. + + private void nearCallback(Object data, DGeom o1, DGeom o2) { + assert (o1 != null); + assert (o2 != null); + if (o1 instanceof DSpace || o2 instanceof DSpace) { + // colliding a space with something + OdeHelper.spaceCollide2(o1, o2, data, nearCallback); + // Note we do not want to test intersections within a space, + // only between spaces. + return; + } + + final int N = 32; + DContactBuffer contacts = new DContactBuffer(N); + int n = OdeHelper.collide(o1, o2, N, contacts.getGeomBuffer()); + if (n > 0) { + for (int i = 0; i < n; i++) { + DContact contact = contacts.get(i); + contact.surface.slip1 = 0.7; + contact.surface.slip2 = 0.7; + contact.surface.mode = dContactSoftERP | dContactSoftCFM | dContactApprox1 | dContactSlip1 | dContactSlip2; + contact.surface.mu = 500.0; // was: dInfinity + contact.surface.soft_erp = 0.50; + contact.surface.soft_cfm = 0.03; + DJoint c = OdeHelper.createContactJoint(world, contactgroup, contact); + c.attach( + contact.geom.g1.getBody(), + contact.geom.g2.getBody() + ); + } + } + } + + + // start simulation - set viewpoint + @Override + public void start() { + //dAllocateODEDataForThread(dAllocateMaskAll); + dsSetViewpoint(xyz, hpr); + System.err.println("Press SPACE to reset the simulation."); + } + private static final float[] xyz = { // [ 3] ={ + -8, 0, 5 + }; + private static final float[] hpr = { // [ 3] ={ + 0.0f, -29.5000f, 0.0000f + }; + + + private static void reset() { + DQuaternion q = new DQuaternion(); + q.setIdentity(); + mbody.setPosition(0, 0, 0 + H); + mbody.setQuaternion(q); + mbody.setLinearVel(0, 0, 0); + mbody.setAngularVel(0, 0, 0); + mbody.enable(); + for (int i = 0; i < halton_numc; ++i) { + DBody body = hbody[i]; + if (body == null) continue; + body.setPosition(halton_pos[i][0], halton_pos[i][1], halton_pos[i][2] + H); + body.setQuaternion(q); + body.setLinearVel(0, 0, 0); + body.setAngularVel(0, 0, 0); + body.enable(); + } + } + + + // called when a key pressed + @Override + public void command(char cmd) { + switch (cmd) { + case ' ': + reset(); + break; + default: + break; + } + } + + @Override + public void step(boolean pause) { + double simstep = 1 / 240.0; + double dt = dsElapsedTime(); + + int nrofsteps = (int) Math.ceil(dt / simstep); + nrofsteps = nrofsteps > 8 ? 8 : nrofsteps; + + for (int i = 0; i < nrofsteps && !pause; i++) { + OdeHelper.spaceCollide(space, 0, nearCallback); + world.quickStep(simstep); + contactgroup.empty(); + } + + dsSetColor(1, 1, 1); + // Draw the convex objects. + for (int i = 0; i < halton_numc; ++i) { + DGeom geom = hgeom[i]; + DBody body = geom.getBody(); + //const dReal *pos = dBodyGetPosition(body); + //const dReal *rot = dBodyGetRotation(body); + DVector3C pos = geom.getPosition(); + DMatrix3C rot = geom.getRotation(); + dsDrawConvex + ( + pos, rot, + halton_planes[i], + halton_numf[i], + halton_verts[i], + halton_numv[i], + halton_faces[i] + ); + } + + if (drawpos) { + dsSetColor(1, 0, 0.2); + dsSetTexture(DS_NONE); + final float l = 0.35f; + for (int i = 0; i < halton_numc; ++i) { + DBody body = hbody[i]; + DVector3C pos2 = body.getPosition(); + float[] pos = pos2.toFloatArray(); + float[] x0 ={ + pos[0] - l, pos[1], pos[2] + } ; + float[] x1 ={ + pos[0] + l, pos[1], pos[2] + } ; + float[] y0 ={ + pos[0], pos[1] - l, pos[2] + } ; + float[] y1 ={ + pos[0], pos[1] + l, pos[2] + } ; + float[] z0 ={ + pos[0], pos[1], pos[2] - l + } ; + float[] z1 ={ + pos[0], pos[1], pos[2] + l + } ; + dsDrawLine(x0, x1); + dsDrawLine(y0, y1); + dsDrawLine(z0, z1); + } + } + } + + public static void main(final String[] args) { + new DemoConvex().demo(args); + } + + private void demo(String[] args) { + DMass m = OdeHelper.createMass(); + + // setup pointers to drawstuff callback functions + // dsFunctions fn; + // fn.version = DS_VERSION; + // fn.start = &start; + // fn.step = &simLoop; + // fn.command = &command; + // fn.stop = 0; + // fn.path_to_textures = DRAWSTUFF_TEXTURE_PATH; + + // create world + OdeHelper.initODE2(0); + world = OdeHelper.createWorld(); + space = OdeHelper.createHashSpace(); + ((DHashSpace)space).setLevels(-3, 5); + OdeHelper.createPlane(space, 0, 0, 1, 0); // Add a ground plane. + + contactgroup = OdeHelper.createJointGroup(); + world.setGravity(0, 0, -9.8); + world.setQuickStepNumIterations(32); + world.setContactMaxCorrectingVel(40); + world.setMaxAngularSpeed(62.8); + world.setERP(0.7); + world.setQuickStepW(0.75); // For increased stability. + + world.setAutoDisableFlag(true); + world.setAutoDisableLinearThreshold(0.01); + world.setAutoDisableAngularThreshold(0.03); + world.setAutoDisableTime(0.15f); + + final float kernelrad = 0.7f; + + mbody = OdeHelper.createBody(world); + mbody.setPosition(0, 0, 0 + H); + m.setSphere(5, kernelrad ); + mbody.setMass( m ); + + for (int i = 0; i < halton_numc; ++i) { + DGeom geom = OdeHelper.createConvex( + space, + halton_planes[i], + halton_numf[i], + halton_verts[i], + halton_numv[i], + halton_faces[i] + ); + hgeom[i] = geom; + final double x = halton_pos[i][0]; + final double y = halton_pos[i][1]; + final double z = halton_pos[i][2]; + final double dsqr = x * x + y * y + z * z; + + if (dsqr < kernelrad * kernelrad && solidkernel) { + geom.setBody(mbody); + geom.setOffsetPosition(x, y, z); + } else { + DBody body = OdeHelper.createBody(world); + hbody[i] = body; + body.setPosition(x, y, z + H); + double volu = halton_volu[i]; + double rad = Math.pow(volu * 3 / (4 * M_PI), (1 / 3.0)); + m.setSphere( 5, rad ); + body.setMass( m ); +//#if 1 + body.setLinearDamping(0.0005); + body.setAngularDamping(0.0300); +//#edif + geom.setBody(body); + } + } + + // run simulation + final int w = 1280; + final int h = 720; + dsSimulationLoop(args, w, h, this); + + contactgroup.empty (); + contactgroup.destroy (); + space.destroy (); + world.destroy (); + OdeHelper.closeODE(); + } + + @Override + public void stop() { + // Nothing + } + +} diff --git a/demo/src/main/java/org/ode4j/demo/DemoFeedback.java b/demo/src/main/java/org/ode4j/demo/DemoFeedback.java index 80af8685..d648e895 100644 --- a/demo/src/main/java/org/ode4j/demo/DemoFeedback.java +++ b/demo/src/main/java/org/ode4j/demo/DemoFeedback.java @@ -89,7 +89,6 @@ private void nearCallback (Object data, DGeom o1, DGeom o2) if ( o1 instanceof DSpace || o2 instanceof DSpace ) { - System.err.println("testing space " + o1 + " " + o2); // colliding a space with something OdeHelper.spaceCollide2(o1,o2,data,nearCallback); // Note we do not want to test intersections within a space, @@ -120,8 +119,8 @@ private void nearCallback (Object data, DGeom o1, DGeom o2) } - private static float[] xyz = { -6, 8, 6}; - private static float[] hpr = { -65.0f, -27.0f, 0.0f}; + private static final float[] xyz = { -6, 8, 6}; + private static final float[] hpr = { -65.0f, -27.0f, 0.0f}; // start simulation - set viewpoint @Override public void start() @@ -160,7 +159,7 @@ private void drawGeom (DGeom g) private static void inspectJoints() { - final double forcelimit = 2000.0; + final double forcelimit = 4000.0; int i; for (i=0; i 1e-6) dDebug(0, "bad u1/u2"); - n.eqCross(u1, u2); - //for (k=0; k<3; k++) tmp[k] = v2[k]-v1[k]; - tmp.eqDiff(v2, v1); - double d = -n.dot(p1); - if (dFabs(n.dot(p1) + d) > 1e-8) dDebug(0, "bad n wrt p1"); - if (dFabs(n.dot(p2) + d) > 1e-8) dDebug(0, "bad n wrt p2"); - if (dFabs(n.dot(p3) + d) > 1e-8) dDebug(0, "bad n wrt p3"); - double alpha = -(d + n.dot(v1)) / n.dot(tmp); - //for (k=0; k<3; k++) tmp[k] = v1[k]+alpha*(v2[k]-v1[k]); + + double error; + // #ifdef dSINGLE + // const dReal uEpsilon = 1e-5, pEpsilon = 1e-6, tmpEpsilon = 1.5e-4; + // #else + final double uEpsilon = 1e-6, pEpsilon = 1e-8, tmpEpsilon = 1e-6; + // #endif + + error = dFabs(dCalcVectorDot3(u1, u2)); + if (error > uEpsilon) dDebug(0, "bad u1/u2"); + + dCalcVectorCross3(n, u1, u2); + + // for (k=0; k < 3; k++) tmp[k] = v2[k] - v1[k]; tmp.eqDiff(v2, v1); - tmp.eqSum(v1, tmp.scale(alpha)); - if (dFabs(n.dot(tmp) + d) > 1e-6) dDebug(0, "bad tmp"); + + double d = -dCalcVectorDot3(n, p1); + + error = dFabs(dCalcVectorDot3(n, p1) + d); + if (error > pEpsilon) dDebug(0, "bad n wrt p1"); + + error = dFabs(dCalcVectorDot3(n, p2) + d); + if (error > pEpsilon) dDebug(0, "bad n wrt p2"); + + error = dFabs(dCalcVectorDot3(n, p3) + d); + if (error > pEpsilon) dDebug(0, "bad n wrt p3"); + + double alpha = -(d + dCalcVectorDot3(n, v1)) / dCalcVectorDot3(n, tmp); + // for (k=0; k < 3; k++) tmp[k] = v1[k] + alpha * (v2[k] - v1[k]); + tmp.eqDiff(v2, v1).scale(alpha).add(v1); + + error = dFabs(dCalcVectorDot3(n, tmp) + d); + if (error > tmpEpsilon) dDebug(0, "bad tmp"); if (alpha < 0) return false; if (alpha > 1) return false; + //for (k=0; k<3; k++) tmp[k] -= p1[k]; tmp.sub(p1); double a1 = u1.dot(tmp); double a2 = u2.dot(tmp); if (a1 < 0 || a2 < 0 || a1 > d1 || a2 > d2) return false; + return true; } diff --git a/demo/src/test/java/org/ode4j/tests/DemoJointsTest.java b/demo/src/test/java/org/ode4j/tests/DemoJointsTest.java index 709924f8..e6afa5b9 100644 --- a/demo/src/test/java/org/ode4j/tests/DemoJointsTest.java +++ b/demo/src/test/java/org/ode4j/tests/DemoJointsTest.java @@ -87,7 +87,7 @@ public class DemoJointsTest { private static DUniversalJoint jointUniversal; - private static int cmd_occasional_error = 0; // perturb occasionally + private static int cmd_occasional_error = 1; // perturb occasionally // info about the current test diff --git a/demo/src/test/java/org/ode4j/tests/DemoOdeTest.java b/demo/src/test/java/org/ode4j/tests/DemoOdeTest.java index 1e7236c3..6402ccd3 100644 --- a/demo/src/test/java/org/ode4j/tests/DemoOdeTest.java +++ b/demo/src/test/java/org/ode4j/tests/DemoOdeTest.java @@ -47,8 +47,28 @@ import static org.ode4j.ode.internal.ErrorHandler.*; import static org.ode4j.ode.internal.ErrorHandler.dSetDebugHandler; - public class DemoOdeTest { + //**************************************************************************** + // matrix sizes + + // #define dALIGN_SIZE(buf_size, alignment) (((buf_size) + (alignment - 1)) & (int)(~(alignment - 1))) // Casting the mask to int ensures sign-extension to larger integer sizes + // #define dALIGN_PTR(buf_ptr, alignment) ((void *)(((duintptr)(buf_ptr) + ((alignment) - 1)) & (int)(~(alignment - 1)))) // Casting the mask to int ensures sign-extension to larger integer sizes + // + // #define MSIZE 21 + // #define MSIZE4 dALIGN_SIZE(MSIZE, 4) // MSIZE rounded up to 4 + private static final int MSIZE = 21; + private static final int MSIZE4 = 24; // MSIZE rounded up to 4 + + + + + + //**************************************************************************** + // matrix accessors + + // #define _A(i,j) A[(i)*4+(j)] + // #define _I(i,j) I[(i)*4+(j)] + // #define _R(i,j) R[(i)*4+(j)] //**************************************************************************** @@ -268,12 +288,6 @@ public void testPlaneSpace() { //**************************************************************************** // test matrix functions - //#define MSIZE 21 - //#define MSIZE4 24 // MSIZE rounded up to 4 - private static final int MSIZE = 21; - private static final int MSIZE4 = 24; // MSIZE rounded up to 4 - - @Test public void testMatrixMultiply() { // A is 2x3, B is 3x4, B2 is B except stored columnwise, C is 2x4 @@ -352,8 +366,10 @@ public void testSmallMatrixMultiply() { @Test public void testCholeskyFactorization() { - double[] A = new double[MSIZE4 * MSIZE], B = new double[MSIZE4 * MSIZE]; - double[] C = new double[MSIZE4 * MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + double[] A = new double[matrixSize]; + double[] B = new double[matrixSize]; + double[] C = new double[matrixSize]; double diff; dMakeRandomMatrix(A, MSIZE, MSIZE, 1.0); dMultiply2(B, A, A, MSIZE, MSIZE, MSIZE); @@ -368,7 +384,8 @@ public void testCholeskyFactorization() { @Test public void testCholeskyFactorizationM3() { - DMatrix3 A = new DMatrix3(), B = new DMatrix3(); + DMatrix3 A = new DMatrix3(); + DMatrix3 B = new DMatrix3(); DMatrix3 C = new DMatrix3(); double diff; dMakeRandomMatrix(A, 1.0); @@ -384,8 +401,10 @@ public void testCholeskyFactorizationM3() { @Test public void testCholeskySolve() { - double[] A = new double[MSIZE4 * MSIZE], L = new double[MSIZE4 * MSIZE]; - double[] b = new double[MSIZE], x = new double[MSIZE], btest = new double[MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + int vectorSize = MSIZE; + double[] A = new double[matrixSize], L = new double[matrixSize]; + double[] b = new double[vectorSize],x = new double[vectorSize],btest = new double[vectorSize]; double diff; // get A,L = PD matrix @@ -445,8 +464,9 @@ public void testCholeskySolveM3() { @Test public void testInvertPDMatrix() { int i, j, ok; - double[] A = new double[MSIZE4 * MSIZE], Ainv = new double[MSIZE4 * MSIZE]; - double[] I = new double[MSIZE4 * MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + double[] A = new double[matrixSize], Ainv = new double[matrixSize]; + double[] I = new double[matrixSize]; dMakeRandomMatrix(A, MSIZE, MSIZE, 1.0); dMultiply2(Ainv, A, A, MSIZE, MSIZE, MSIZE); @@ -501,7 +521,8 @@ public void testInvertPDMatrixM3() { @Test public void testIsPositiveDefinite() { - double[] A = new double[MSIZE4 * MSIZE], B = new double[MSIZE4 * MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + double[] A = new double[matrixSize], B = new double[matrixSize]; dMakeRandomMatrix(A, MSIZE, MSIZE, 1.0); dMultiply2(B, A, A, MSIZE, MSIZE, MSIZE); assertFalse(dIsPositiveDefinite(A, MSIZE)); @@ -520,9 +541,11 @@ public void testIsPositiveDefiniteM3() { @Test public void testFastLDLTFactorization() { int i, j; - double[] A = new double[MSIZE4 * MSIZE], L = new double[MSIZE4 * MSIZE]; - double[] DL = new double[MSIZE4 * MSIZE]; - double[] ATEST = new double[MSIZE4 * MSIZE], d = new double[MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + int vectorSize = MSIZE; + double[] A = new double[matrixSize], L = new double[matrixSize]; + double[] DL = new double[matrixSize]; + double[] ATEST = new double[matrixSize], d = new double[vectorSize]; double diff; dMakeRandomMatrix(A, MSIZE, MSIZE, 1.0); dMultiply2(L, A, A, MSIZE, MSIZE, MSIZE); @@ -546,9 +569,11 @@ public void testFastLDLTFactorization() { @Test public void testSolveLDLT() { - double[] A = new double[MSIZE4 * MSIZE], L = new double[MSIZE4 * MSIZE]; - double[] d = new double[MSIZE], x = new double[MSIZE]; - double[] b = new double[MSIZE], btest = new double[MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + int vectorSize = MSIZE; + double[] A = new double[matrixSize], L = new double[matrixSize]; + double[] d = new double[vectorSize], x = new double[vectorSize]; + double[] b = new double[vectorSize], btest = new double[vectorSize]; double diff; dMakeRandomMatrix(A, MSIZE, MSIZE, 1.0); @@ -571,9 +596,11 @@ public void testSolveLDLT() { @Test public void testLDLTAddTL() { int i, j; - double[] A = new double[MSIZE4 * MSIZE], L = new double[MSIZE4 * MSIZE]; - double[] d = new double[MSIZE], a = new double[MSIZE]; - double[] DL = new double[MSIZE4 * MSIZE], ATEST = new double[MSIZE4 * MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + int vectorSize = MSIZE; + double[] A = new double[matrixSize], L = new double[matrixSize]; + double[] d = new double[vectorSize], a = new double[vectorSize]; + double[] DL = new double[matrixSize], ATEST = new double[matrixSize]; double diff; dMakeRandomMatrix(A, MSIZE, MSIZE, 1.0); @@ -609,15 +636,17 @@ public void testLDLTAddTL() { @Test public void testLDLTRemove() { int i, j, r; - int[] p = new int[MSIZE]; - double[] A = new double[MSIZE4 * MSIZE]; - double[] L = new double[MSIZE4 * MSIZE]; - double[] d = new double[MSIZE]; - double[] L2 = new double[MSIZE4 * MSIZE]; - double[] d2 = new double[MSIZE]; - double[] DL2 = new double[MSIZE4 * MSIZE]; - double[] Atest1 = new double[MSIZE4 * MSIZE]; - double[] Atest2 = new double[MSIZE4 * MSIZE]; + int matrixSize = MSIZE4 * MSIZE; + int vectorSize = MSIZE; + int[] p = new int[vectorSize]; + double[] A = new double[matrixSize]; + double[] L = new double[matrixSize]; + double[] L2 = new double[matrixSize]; + double[] DL2 = new double[matrixSize]; + double[] Atest1 = new double[matrixSize]; + double[] Atest2 = new double [matrixSize]; + double[] d = new double[vectorSize]; + double[] d2 = new double[vectorSize]; double diff, maxdiff; // make array of A row pointers diff --git a/demo/src/test/java/org/ode4j/tests/TestIssue0008_Gimpact.java b/demo/src/test/java/org/ode4j/tests/TestIssue0008_Gimpact.java index 7cf8604f..3a0077a5 100644 --- a/demo/src/test/java/org/ode4j/tests/TestIssue0008_Gimpact.java +++ b/demo/src/test/java/org/ode4j/tests/TestIssue0008_Gimpact.java @@ -57,8 +57,7 @@ public class TestIssue0008_Gimpact { */ @Test public void testLargeTrimesh() { - //TODO Test does not fail yet, even if the last part (0x3FF) is removed from line 76 in - //RadixSort. + //TODO Test does not fail yet, even if the last part (0x3FF) is removed from line 76 in RadixSort. // float[] size = new float[]{ 5.0f, 5.0f, 2.5f }; //